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ABSTRACT 


Wireless communication technology has become a critical aspect in many civilian and 
military applications. With regard to remote sensing, search and rescue, disaster relief 
operations and signals intelligence, there exists an interest in developing capabilities to 
collect these signals-of-interest. The objective of this dissertation is to maximize signal 
collection performance in the presence of signal measurement and sensor related errors. 
To accomplish this objective, we proposed a signal collection scheme that exploits an 
elevated, mobile network to maximize the collaborative collection of a target signal. 

The proposed scheme begins with source localization. This technique consists of 
an initial weighted least-squares estimate followed by a maximum-likelihood estimate. 
Implemented on an elevated, mobile network, this technique is able to obtain an optimal 
localization. To enhance localization robustness, we developed an outlier rejection 
process that mitigates the effects of measurement and sensor position errors. 

To collect the signal, this research quantified the effects of sensor position errors 
on beamforming and proposed a novel signal collection scheme that combines signal 
estimation and collaborative beamforming. Using all these techniques in concert, we 
were able to show that the proposed scheme outperforms standard collaborative 


beamforming in the presence of sensor position errors. 
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LIST OF ACRONYMS AND ABBREVIATIONS 


ASLA airborne symmetric line array 
AWGN additive white Gaussian noise 
BPSK binary phase-shift keying 
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RDOA range difference-of-arrival 
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EXECUTIVE SUMMARY 


Wireless communication technology has penetrated many aspects of both civilian and 
military applications. With cellular technology constantly evolving and proliferating, we 
now see wireless communications as an integral part of daily life. With regard to remote 
sensing, search and rescue, disaster relief operations and signals intelligence, there exists 
an interest in developing capabilities to collect these signals. The objective of this 
dissertation is to maximize signal collection performance in the presence of various 
signal- and sensor-related errors. To accomplish this objective, we proposed a signal 
collection scheme that exploits an elevated, mobile network to maximize the 


collaborative collection of a target signal. 


With its roots dating back to the development of radar during World War II, 
beamforming has found a prominent place in modern wireless communications. In its 
classical use, through the manipulation of array weights, beamforming coherently 
amplifies a signal in a given direction, while reducing undesirable signals in others. 
Although, beamforming is a powerful technique that enables spatial filtering, it has a 
strong dependence on the a priori knowledge of the target signal’s direction of arrival 


[1]. Because of this dependency, beamforming is often initialized by source localization. 


Our proposed scheme begins with a source localization technique; this technique 
is used to determine the signal’s position and direction of arrival. This knowledge is then 
used by a collaborative beamformer to maximize and collect the target signal. Both of 
these techniques are enhanced through the use of an elevated, mobile network. This 
mobility allows for the reconfiguration of the sensor network’s topology to create an 
ideal sensor-target geometry. This geometry minimizes any geometric dilutions of 
precision, thus allowing for an optimal location estimate to be achieved [2]. The 
network’s elevation increases signal power and range while enabling unique sensor 


formations for signal collection that is robust against sensor position errors. 


The proposed scheme begins with a localization technique. In this localization, 


we propose a two-stage technique capable of obtaining an accurate localization in the 


XVil 


presence of signal noise and sensor position errors. This technique consists of an initial 
weighted least-squares estimate followed by a refining maximum-likelihood estimate. 
Here, the initial estimate is used to reconfigure the network’s topology. This 
reconfiguration is used to create an optimal sensor network to target geometry, thus 
creating the necessary conditions for the refining estimation to deliver an optimal location 
estimate. Our simulation results show the proposed localization technique to be efficient, 


i.e., 1t approaches the Cramer-Rao lower bound in the small error region [2]. 


To enhance the localization performance, we developed a measurement outlier 
rejection process to mitigate the effects of measurement and sensor position errors. This 
technique uses a combination of single case diagnostics and the Mahalanobis distance to 
identify and remove specious measurements for the least-squares estimate [3]. Through 
simulation, we demonstrated this technique to be effective in the presence of both 


measurement and position errors. 


To develop the proposed signal collection scheme, we began by analyzing the 
effects of sensor position errors on the array factor. We derived an expression for the 
effects of position error on the array’s main beam gain. Through simulation, we 
validated our results and showed that the mean value of the main beam signal phase was 
unaffected by position errors. To enhance the collection performance in the presence of 
these errors, we developed a signal estimator for both uniform and Gaussian position 
errors. Our simulation results showed a maximum improvement in array gain of 
approximately 37 percent for the standard deviation of position error values greater than 
0.4 m. Other results demonstrated that our approach of combining a signal estimation 
technique with collaborative beamforming is a viable means of collecting a target signal 


from an elevated, mobile wireless sensor network. 


Taking advantage of the concept of a unique elevated, mobile wireless sensor 
network realized through the use of multirotor UAVs, our scheme used two existing 
localization techniques to deliver a precise source location estimate. Using this 
information with statistical knowledge of the sensor position errors in a signal estimator, 
we were able to derive a novel collaborative signal collection scheme. This scheme was 


shown to be capable of collecting and amplifying a target signal in the presence of such 
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errors. With all these techniques in concert, the objective of this dissertation to maximize 
signal collection performance in the presence of various signal- and sensor-related errors 


was achieved. 
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I. INTRODUCTION 


Wireless communication technology has penetrated many aspects of both civilian 
and military applications. The International Telecommunications Union (ITU) estimates 
that there are approximately seven billion cellular subscribers worldwide [1]. This 
translates to a global penetration rate of 95.5 percent and is expected to increase each 
year. Clearly, wireless communications has become a common aspect of human life. 
With regard to remote sensing, search and rescue, disaster relief operations and/or signals 
intelligence, there exists an interest in developing capabilities to collect these signals. 


One such solution is collaborative beamforming [2]. 


With its roots dating back to the development of radar during World War II, 
beamforming has found a prominent place in modern wireless communications. 
Beamforming is much akin to the matched filter but from a spatial perspective [3]. 
Instead of linearly weighting a time domain signal for signal filtering, a beamformer 
linearly weights spatial array data for spatial filtering. In its classical use, beamforming 
coherently amplifies a signal in a given direction while reducing undesirable signals in 
other directions. Although beamforming is a powerful technique that enables spatial 
filtering, it has a strong dependence on a priori knowledge of the target signal’s direction 
of arrival [4]. Because of this sensitivity, beamforming is often preceded by a direction 


of arrival estimation phase. 


With regard to direction-of-arrival estimates, hyperbolic localization [5] has 
appealing qualities of interest to this research. Of particular interest is hyperbolic 
localization performed from an elevated, mobile wireless sensor network (WSN). The 
term hyperbolic is used to denote the use of signal time-difference-of-arrival 
measurement and calculations. The network is elevated to promote a line-of-sight 
communication link for better range and reliability, while the mobility allows for the 
reconfiguration of sensor-target geometries for more precise direction-of-arrival 


estimates. 


The extension of beamforming to wireless sensor networks has given rise to the 
concept of collaborative beamforming. Here, a fully synchronized network of sensors is 
used to coherently collect and combine multiple received signals. These signals along 
with network topology information are used to form a distributed beamforming array. 
This technique coupled with recent advances in micro sensor and multirotor unmanned 
aerial vehicle (UAV) technologies has led to the concept of a semi-stationary airborne 
sensor network. Such a network is capable of providing both the elevation and mobility 
needed for highly accurate hyperbolic localization while providing a suitable network 


platform for collaborative signal collection. 


A. OBJECTIVE 


The dissertation objective is to maximize signal collection performance in the 
presence of various signal and sensor related errors. To accomplish this task, we propose 


a new scheme for signal collection that is preceded by a localization phase. 


To improve the scheme’s localization accuracy, we propose the use of an 
optimally positioned two-stage localization technique. This technique consists of an 
initial weighted least-squares estimate followed by a refined maximum-likelihood 
estimate. Here, the initial estimate is used to reconfigure the network’s topology, thereby 
creating an optimally configured network to obtain the refined hyperbolic estimate. To 
combat the effects of noise, we propose the use of measurement outlier diagnostics in 
concert with the two-stage technique. Here, we postulate that the identification and 
removal of specious time-difference-of-arrival measurements will improve the robustness 


of localization in the presence of noise. 


To maximize the signal collection performance in the presence of noise, we 
propose the use of signal estimation in conjunction with collaborative beamforming. 
Using the concept of sensor stacks with knowledge of the sensor positional error 
statistics, we can derive a signal estimator to combat the effects of sensor position errors. 
To further enhance the collection capability for noise and interference rejection, we also 


propose the use of virtual filling to increase the array’s ability to reject interfering signals. 


B. RELATED WORK 


Beamforming has been shown to be an effective method for signal collection and 
interference rejection [6], but it has been shown to be highly susceptible to array steering 
vector errors [7]. This is when an array steering vector is not in line with that of the 
target signal. To make beamforming more robust against these array mismatch errors, 
Ahmed and Evans [8] suggest the use of inequality constraints on the array weights while 
Er and Cantoni [9] suggest the use of derivative constraints. Both of these approaches, 
although effective, also decrease the array’s interference rejection capabilities. Li, Stoica 
and Wang [10] developed a robust beamformer that minimized the effects of non-random 
steering errors using a diagonal loading method. Lee and Lee [11] proposed a robust 
beamformer for signal collection, which minimizes a cost function based on received 
signal data and knowledge of steering error statistics. The signal collection scheme 
proposed in this research uses elevated, mobile sensor network formations and signal 
estimation in the context of collaborative beamforming. Preceded by a localization phase 
that creates a unique sensor network-to-target formation, the scheme is able to isolate the 
effects of phase perturbations due to position errors. It then compensates for them 
through the use of sensor formations, signal estimation, and knowledge of sensor position 
error statistics. Unlike the schemes of Ahmed and Evans [8] and Er and Cantoni [9], our 
scheme does not require any constraints on the beamformer weights. Similar to that of 
Lee and Lee [11], our scheme utilizes sampled array data and knowledge of error 
statistics to compensate for array phase errors. The difference being that our scheme uses 
the sampled data in a signal estimation context with consideration to array formation 
rather than an array weight optimization problem with no consideration to array 


formation. 


One of the earliest hyperbolic location systems was investigated by Carson [12] in 
1972. He proposed a system of hydrophones to collect signal arrival times, with a goal to 
estimate the location of subsurface acoustic emitters. Today, numerous extensions and 
breakthroughs have been made to adapt these concepts for radio frequency (RF) signals. 
Fang [13] provides a closed form solution, but the estimator is not efficient, i.e., it does 


not achieve the Cramer-Rao lower bound [14], [15]. Hahn and Tretter [16] used the 
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hyperbolic equations in a maximum-likelihood estimation scheme that is unbiased and 
asymptotically efficient. Being a linearized system of non-linear equations, the 
maximum-likelihood estimator requires an initial reference estimate to avoid local 
minima. Chan and Ho [5] propose a weighted least-squares approach that yields an 
efficient estimator in the small error region. Our approach will combine the maximum- 
likelihood estimate of Hahn and Tretter [16] and the weighted least-squares estimate of 
Chan and Ho [5] to create a two-stage localization technique. Furthermore, this 
technique will use optimal sensor network formations as outlined by Ho and Vicente 
[17]. This combination of hyperbolic estimation and optimal sensor network formations 
will be used to provide highly accurate localization from an elevated, mobile WSN 


platform. 


Rousseeuw and Leroy [18] provided the first comprehensive review of robust 
estimation and the effects of statistical outlier measurements. Key contributions by 
Mahalanobis [19] and Cook [20] helped quantify the effects of outliers in estimation by 
defining the Mahalanobis and Cook distance. Picard and Weiss [21] coupled these 
findings with the concept of sparse representation to derive a bound on the number of 
outliers that can be identified in a hyperbolic localization problem. Many researchers 
have contributed techniques to combat the effects of outlier measurements. McGuire, 
Plataniotis and Venetsanopoulos [22] proposed the use of data fusion techniques. Yang, 
Wang and Luo [23] showed that convex relaxation methods can be used to overcome 
outliers induced by both measurement noise and sensor position errors. Choi, et al. [24] 
used recursive filtering to deal with errors from incorrect stochastic information in 
location estimators. With the proposed localization technique implemented on an 
elevated, mobile sensor network, there will be increased measurement noise due to sensor 
position errors. To minimize the effects of these errors as well as measurement outliers, 
we propose an outlier rejection process based on single case diagnostics [18] and the 
squared-Mahalanobis distance [19]. Overall, a robust localization scheme that combines 
the proposed two-stage localization technique and the outlier rejection process is 


achieved. 


C; ORGANIZATION 
This dissertation is organized as follows. 


The fundamental principles involved with hyperbolic localization and 
beamforming is the focus of Chapter II. Localization based on the time-difference-of- 
arrival measurements, weighted least-squares estimator, maximum-likelihood estimator, 
and robust signal outlier detection techniques are discussed in the localization section. 
With collaborative signal collection as the end goal, we also provide an introduction to 


beamforming and its extension to collaborative beamforming. 


The proposed two-phase scheme is presented in Chapter HI. This includes a 
description of the scheme’s process flow as it moves through the two phases of 
localization and collection. A discussion of the deployment and operating concepts is 
provided to give context to the proposed scheme. The many figures of merit used to 


evaluate the scheme’s performance are also detailed. 


The detailed analysis in support of the robust localization phase of our scheme is 
included in Chapter IV. Additional robust signal processing techniques that are 
implemented in both estimators to combat the effects of measurement outliers and sensor 
position errors are also included. The performance of the robust localization phase is 


then examined using simulation. 


The signal collection phase of the proposed scheme is described in Chapter V. 
This includes an analysis of sensor position errors and their effects on an array’s array 
factor. Based on these findings, we derive a novel signal estimator to combat the effects 
of uniform and Gaussian sensor position errors. The performance of the resulting 


collection scheme is then investigated using simulation. 


In conclusion, a summary of the dissertation, a discussion of the research’s 
significant contributions, and a list of potential future research topics are given in Chapter 


VI. An appendix is also included, which contains selected MATLAB simulation scripts. 
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Il. BACKGROUND 


To perform signal collection, we must understand the different processes involved 
in source localization and collection. In this section, we provide the necessary 
background material to properly present our research and results. We begin with a brief 
discussion on WSN that forms the basis of the proposed signal collection scheme. With 
the premise of sensor networks explored, we discuss the principles of hyperbolic/TDOA 


source localization and collaborative beamforming. 


Throughout the dissertation, the following notation is used. Bold lower-case font 


indicates a vector, while bold upper-case font indicates a matrix, Pr(*) indicates the 


a Ws Tks as 
probability of an event, £ {x} indicates the expected value, (*) indicates the transpose 





operation, (*) indicates the Hermitian transpose operation, | indicates the 2-norm 


operation, I is the identity matrix, 1; indicates a column vector of length i, tr(*) 
indicates the trace of a matrix, j indicates an imaginary quantity, c is the signal 


propagation speed, and 1 is the signal wavelength. 


A. WIRELESS SENSOR NETWORKS 


Advances in wireless communications coupled with low cost multifunctional 
sensors have given rise to the growth of WSNs. A WSN is a network comprised of many 
sensor nodes deployed in or near a phenomenon to be observed. Typical 
implementations have identical sensor nodes collecting data and routing it through a 
shared channel to a sink node for processing/storage. An illustration of this concept is 
shown in Figure 1. Often, the deployment mechanisms are not predetermined, thus 
requiring arbitrary deployment schemes in inaccessible locations. The sensor network 
nodes collaborate to achieve a common goal. Some examples of this collaboration can be 
found in [25]. For example, sensor nodes can be distributed across a geographic region 
to monitor rainfall to support flood predictions, covert sensors could be deployed along 


known supply lines to monitor enemy troop movements, or RFID sensors could be 


deployed throughout a warehouse for inventory control [25]. A few other notable 
research areas where WSNs have been applied are environmental monitoring [26], 
military applications such as acoustic detection of helicopters and chemical weapons 


detection [27], and logistical applications such as global shipment tracking [28]. 





Sensor 
Field 





Figure 1. Wireless sensor network concept. 


1. Ground-based, Stationary Wireless Sensor Network System 
Architecture 


Sensor networks have many benefits, but the realization of an actual network can 
be complex with many design factors to consider. Some key factors are fault tolerance, 
scalability, production cost, hardware constraints, sensor network topology, operating 


environment, transmission media, power consumption, and protocol stack [29]. 


A major aspect of WSNs is the hardware components of each sensor node. A 
typical sensor node has four basic components: a sensor, processor, transceiver, and 
power source. There can also be some application dependent components, such as 
localization, power, and mobility systems. Each system plays a key role in this body of 


research. The key components of such a sensor node are shown in Figure 2. 


Processor 


Communication 


Transceiver Storage 


Power Unit 





Figure 2. Key components of a sensor node. 


2. Network Time Synchronization 


Clock synchronization is an underling requirement for many WSN applications. 
A key support function, time synchronization is vital to such applications as data fusion, 
transmission time synchronization, time-based channel scheduling, and node sleep/wake 
scheduling [30]. Most studies on clock synchronization have revolved around protocol 
design [31], but at its core clock synchronization is a parameter estimation problem, thus 


meriting a signal processing approach. The clock in each WSN node can be defined as 


c(t)=t, (1) 


where ¢ is the ideal reference time. Realistic systems deviate from the ideal clock due to 


imperfect clock oscillators, resulting in the general clock function 


c(t) =st+0,, (2) 


where s, is the clock skew and o, is the offset. From (2), we can see that if not 


addressed, clock error can grow linearly over time. For example, a typical crystal-quartz 
oscillator commonly used in WSNs can have frequencies that vary up to 40 ppm, causing 


potential clock errors up to + 40 ws per second [30]. Although, clock synchronization is 


a concern in many WSN applications, this research does not address this problem and 


assumes perfect synchronization across the network. 


3. Elevated, Mobile Wireless Sensor Networks 


It is well known that placing a received antenna on an elevated platform improves 
signal reception. By elevating the antenna, the communication link is less likely to be 
obstructed while promoting a line-of-sight signal path. From [32], the received power at 


an antenna over a reflective surface can be expressed as 





hh | GG 
=| r -P, (3) 
al 





where P is the transmitted power, /, is the receiver antenna height, , is the transmitter’s 
antenna height, d is the distance, G.is the gain of the receive antenna, G,is the gain of 
the transmit antenna, and L, is other associated losses. From this expression, we can see 
that an increase in receiver or transmitter height significantly decreases the signal path 
loss, i.e., increases the ratio of P toP . For example, given a scenario in which 
h.=h=2m, G.=G,=L,=04dB and d=5000m, the resulting path loss is 136 dB. 


By increasing both h, and h, to 100 m, the path loss is decreased to 68 dB. 


Mobility in sensor networks has many advantages with regard to source 
localization and collection. In this research, it is not the velocity of a node that is 
important but its ability to reposition and maintain position. For example, in 
collaborative beamforming, a mobile sensor network can be reconfigured to achieve a 
desirable array factor [2], [33]. In hyperbolic localization, a sensor network can be 
reconfigured to obtain optimal sensor to target formations for increased location 


estimation accuracy [17], [34]. 


To leverage the advantage of both elevation and mobility, the proposed scheme is 


developed in the context of an elevated, mobile WSN; elevated to promote a LOS signal 
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path for increased received power and range, and mobile for control over sensor-target 


geometry for increased localization accuracy. 


4, Multirotor UAV as an Elevated, Nobile WSN Sensor Node 


Sensor nodes for WSNs have been developed for different operating 
environments, such as stationary ground networks [35], ground mobile networks [36], 
and airborne networks [37]. For signal collection, the ideal operating environment is 
elevated above physical obstructions and noisy surface environments. Such a vantage 
point extends communications range and promotes a line-of-sight signal path. For this 
research, we propose the use of an elevated, mobile WSN. A possible realization of this 
network is via the use of multirotor unmanned aerial vehicles (UAVs) [38], [39]. These 
UAVs can provide a mobile sensor platform that in contrast to fixed-wing UAVs can 
maintain a given position. This research does not undertake sensor network 
implementation using multirotor UAVs; however, the simulations in this dissertation 


utilize the concept of an elevated and mobile sensor node realized by multirotor UAVs. 


Although the modern multirotor UAV has only recently been developed, it 
remains a relatively simple machine. For instance, a quadrotor has four rotors held 
together with a rigid frame [40]. Control of such vehicles is basically accomplished 
through the differential manipulation of each rotor’s thrust [38], [39]. In contrast to 
fixed-wing UAVs, multirotor UAVs are under-actuated, with the remaining degrees of 
freedom controlled through system dynamics. Having demonstrated their effectiveness 
in many applications, such as environment mapping and monitoring [41], [42], 
transportation and construction [43], and wireless communication/networking [44], they 
have become a popular research platform. For more information regarding multirotor 


UAVs, see [38], [39]. 


a. Multirotor UAV Station Keeping 


There are many aspects to multirotor UAVs to consider when used as a sensor 
node. Of key concern is its ability to station keep a given position, for any deviation 
from its given position due to system errors and/or wind results in sensor position errors. 


The station keeping performance of a standard proportional-integral-derivative (PID) 
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controlled quadrotor in the presence of wind is found in Figure 3. Here, the wind is 
modeled using a Dryden model [45]. From the illustration, we can see the observed 
position errors due to station keeping operations do not easily conform to any standard 
distribution. The actual position errors of a multirotor UAV are mainly governed by 
wind, GPS accuracy, and flight controller scheme [46], [47], [48]. In the absence of such 


knowledge, we resort to two basic position error distributions, uniform and Gaussian. 
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Figure 3. Simulated quadrotor station keeping in the presence of Dryden modeled 
wind, a) UAV three-dimensional (3D) motion, b) histogram of two- 
dimensional (2D) position errors, from [48]. 


B. HYPERBOLIC SOURCE LOCALIZATION 


Hyperbolic localization is the process of determining the location of a signal 
emitter based on difference of a signal’s time-of-arrival (TOA) across an array of sensor 
nodes. Because of its dependence on the difference in TOA, it is also called localization 
via time-difference-of-arrival (TDOA) measurements. Hyperbolic localization has been 
studied extensively and can be found in many applications [49], [50], [51]. Hyperbolic 
location is commonly applied in WSN as a passive technique; passive in the sense that 
the receiving sensors require little information or synchronization with the transmitter. 


Two-dimensional (2D) hyperbolic localization begins when a signal is received by three 
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or more synchronized sensors [5]. The transmitter’s location is then determined by a set 


of hyperbolic equations that are derived from the TDOA measurements. 


To derive this system of hyperbolic equations, let us consider a signal’s TOA at 


the i" node of an M node WSN as 

behedVe, (4) 
where fo, di, and c are the signal’s time-of-emission, distance from the transmitter to the 
i" node, and signal propagation speed, respectively. Typically, these TOA measurements 


are obtained using the cross-correlation method [16]. By taking the difference between a 


pair of nodes’ measured TOA, we eliminate ¢) and can obtain a set of TDOA equations as 
tia —t, =(4,4,-4,) 7c, (5) 


with i=2,3,...,é@. This can be rewritten to generate hyperbolic equations that include 


the transmitter’s and the i" sensor’s location in Cartesian coordinates as [52] 


e(t.1-t) = (ia 2): + (Yet y,) V(~% 2%.) H(y, y,) (6) 


where (x.,y,) is the 7” node’s location and (x, is the transmitter’s location. With a 
Jj pyr 





minimum of three sensor nodes yielding two simultaneous equations in the form of (6), a 


two-dimensional location estimate can be determined [5]. 


Since the system of equations governing this reconstruction is non-linear, solving 
them becomes a non-trivial matter. There are two ways to deal with this issue. One is to 
avoid the non-linear equations altogether by selectively setting the coordinate system 
[34], [5], and the other linearizes the function through Taylor series expansion about a 
reference point [17], [52]. Additionally, the research in this dissertation is focused on 
elevated, mobile networks with an emphasis on a line-of-sight (LOS) signal path and 


does not address the non-LOS case. For non-LOS solutions, see [53] and [54]. 
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1. Weighted Least-Squares Estimation 


The least-squares method is the standard approach for parameter estimation, given 


a set of measurements that are derived from a linear function of inputs and parameters 
[55]. The least-squares solution is obtained as a set of NV, parameters that minimizes the 


sum of the squared errors: 
Eis = (3, —f (%,@5)) > (7) 


where N,, is the number of noisy measurements, f (x,,@,5) is a linear function of the i” 


input variable x, and the N,, x1 solution parameter vector @,,, and y, is the i” noisy 


measurement defined as 
Ji = f (%,@5)+é,, (8) 


with €, being the zero mean Gaussian measurement error. It is important to clarify the 
linear nature of f(x,,@,,). This function is considered a linear model, when the 


relationship between the input variable and the output is a linear parameterization of @,,. 


When the measurements are taken at different levels of uncertainty, the weighted least- 
squares approach has been shown to be optimal [56], [57]. The weighted least-squares 


solution minimizes the sum of weighted squared-errors: 
Nn 2 
Ems = We (y, — f (%;s@z5)) > (9) 
i=l 


where w,,, is the i" weight and @,,, is the N,x1 parameter vector. The weighted 


least-squares estimate in matrix form is expressed as [56] 
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Oyys =(X™WX)' X™Wy, (10) 


where yis N,, x1 vector of noisy measurements, the J, x N,, weighting matrix is 





— 0 
comm 
We=| : “< : : (11) 
1 
0 2 
L Oy.Nn Al 








where o, is the measurement variance and the NV, x N,,, input variable matrix is [56] 








Of (x,,@,5) ee Of (x,,@,5) 
oo, Oy, 
x= Fs) | : . (12) 
00, 5 
of (sey , 0,5) of (x5 , 0,5) 
[ oo, OWy | 
2. Maximum-Likelihood Estimation 


The maximum-likelihood estimator also requires the system of equations to be 
linear. The term “maximum-likelihood” is used because the solution maximizes the 
likelihood, i.e., the statistical model of the estimate matches that of the measurements. 
Since the maximum-likelihood estimator is both asymptotically unbiased and efficient, 
i.e., achieves the Cramer-Rao lower bound (CRLB) [15], [58], it has become widely 


adopted in the field of parameter estimation. 


To derive the maximum-likelihood solution, we consider an Mx1 noisy 


measurement vector given by 


y=2(@,x)+8, (13) 
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where @ is the N.,x1 vector of unknown but nonrandom set of parameters to be 
estimated, Z (a, x) is a function of @ and the input vector x, and ¢€ is the zero mean 


Gaussian measurement error. The likelihood of y for a given @ is governed by its 


conditional probability density function and is expressed as [52] 


fiw (¥|@) = 


caren aLy elo] Ce[y-e(oy]}, 


é 








where F denotes the determinant of a matrix, and C, is the covariance matrix of the 


measurement error ¢ and is defined as 
C. =E{(Ga-E{6.0})( Gu -F (603) }. (15) 


Finally, the maximum-likelihood estimator of @ can be obtained by minimizing the 


exponent quadratic of (14) [52] 








2. =[y-2(o.x)] C;"[y-z(o.x)]. 


We now encounter the same problem as with the weighted least-squares estimator, where 
z(@, x) is nonlinear, i.e., the relationship between y and x is not a linear parameterization 
of w. The standard solution is then to linearize the functions through a Taylor series 
expansion about a reference point @.. Using only the first two terms of the expansion, 


we have the following approximation [52] 


z,(@,x) *z(@,x)+G,(@-@,), (16) 
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where G, is the MxN,,, gradient matrix given by 














il Zi 
0a, _ Orn, , 7 
G,=| ; (17) 
Hse: ae. Oy, 7 | 


The maximum-likelihood estimator can then be obtained as [52] 


Oy, =0,+(GICIG,) GLC! (y-z(w,,x)). (18) 


ts é 


It is important to note that the maximum-likelihood estimator can be iteratively 
processed if there are multiple sets of y,,. These iterations can be used to overcome the 
effects of an inaccurate @, or compensate for increased measurement noise. The 
iterative performance of a 15 node maximum-likelihood estimator initialized with an 
inaccurate @, is shown in Figure 4. With @. set 1100 meters away from the true 
location, we can see that the estimator quickly approaches the asymptotic optimal 


minimum error value [59]. 
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Figure 4. Iterative performance of a hyperbolic maximum-likelihood estimator in 
the presence of an inaccurate @. estimate. 
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3. Optimal Formations for Source Localization 


To derive the optimal formations for hyperbolic localization, we first determine 
its Cramer-Rao lower bound (CRLB) [17]. For localization using polar coordinates, the 
CRLB is derived using the Fisher information matrix of the hyperbolic localization. The 


Fisher information matrix is given by [17] 
yo  ¥ 0 470 72 0 70 
Vs. (25 or 00 or 06 


where r is range, @is bearing, d°is the (M-1)x1 true TDOA vector, and the 














(M —1)x(M -1) covariance matrix of the TDOA measurement error C, is expressed as 
[17] 


2 


o 
C,= ral +1y My); (20) 
where oy, is the range measurement error variance, I is the (M —1)x(M —1) identity 


matrix, and 1,,, is a column vector of (M —1) ones. The elements of the Fisher 


information matrix from (19) are [17] 


2 Degli A 
2Mo,c'r, ial ial 


Ae sara | ME (0-8.)-| Sirs) ; (21) 


=I ae oe se | 
We = cae MSs (@.-4)-| Sis" (8) [Srsin(—a) |. 


RE I, i=l i i 








and 
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M M : 
2. = aM Sa (a,-0)-[Z in) . (23) 


i=] 


where c is the signal propagation speed and (7,6, ) is the i” sensor’s polar location. 
Based on the use of the maximum-likelihood estimator, which is asymptotically efficient, 
the objective is to choose a sensor formation that will decouple the range and bearing 
estimates, i.e., a formation that makes the element Y,=0. In [17], it is shown that a 
system of concentric circles satisfies this condition. From within this concentric circle 
formation, we can select a sub-formation that either maximizes X, or Z,, 1.¢., 


minimizes the CRLB of range or bearing 


(24) 





In [52], Torrieri showed that the formation of the sensors in relation to the signal 
emitter affects the localization accuracy. This degradation in accuracy is termed the 
geometric dilution of precision (GDOP). For the maximum-likelihood estimator, it is 


given by [52] 





og - |t((ac"6,)') /(co,), (25) 
where 
o = ye (26) 
S Ny = ti 
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is the average variance of the arrival times with o/, being the variance of the TOA 


measurements at the i" node. To combat this degradation, we employ optimal sensor 
formations that minimize the effects of GDOP [17]. Through inspection of the CRLB, it 


can be shown that the sensor formations from Figure 5 minimize GDOP [34], [17]. 

Each of the three formations in Figure 5 is suited to different applications 
depending on the estimated parameter of interest. Here, J is the number of sensor groups 
and M is the total number of sensors. The circular formation(J =3orJ >5) provides 
optimal bearing and range estimates, independent of the targets bearing. The line 


formation (J = 2) delivers the best range and bearing estimates but is dependent on the 


target’s bearing [17], [34]. The cross formation (J = 4) delivers performance comparable 


to the circular formation with minor dependence on the target’s bearing. Although some 
of the formations in Figure 5 may not look like concentric circles, like the linear or cross 
formations, they are indeed two concentric circles. One large circle surrounding a 


smaller one whose radius is zero. 


Linear (J = 2) Cross (J = 4) 
M/3 


Circular 


M/s (J=3 or J25) 





Figure 5. Optimum array geometries for position estimation (minimum 
uncertainty area) using TDOA localization. J = number of sensor groups, @ 
= total number of sensor, from [17]. 


For more information regarding these optimal formations see [17] and [34]. 
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4. The Airborne Symmetrical Line Array Network Configuration 


In this section, we describe the term airborne symmetric line array (ASLA) 
formation. The ASLA formation consists of multiple elevated, mobile sensors in a line 
formation symmetric about the origin. This formation is based on the linear formation in 
Figure 5 and is illustrated in Figure 6. It consists of three sensor stacks that contain 


N,=M /3 sensors each. The center stack is located at the origin, with two outer stacks 
flanking the center at a distance of r,. For hyperbolic localization, this formation has 


been show to minimize the maximum-likelihood error variance [17]. Note that each 
sensor in a stack must be separated by 4/2 to guarantee independent noise and to 


minimize mutual coupling [34]. 


X-axis 
‘ Outer stack Center stack Outer stack 


Figure 6. Illustration of the ASLA formation. 


The linear formation is special in the sense that it delivers the best performance in 
bearing and range estimates. The caveat is that this optimality is dependent on the target 
emitter’s bearing [34]. A Monte Carlo simulation illustrating the dependence on the 
target bearing is shown in Figure 7. Here, the hyperbolic localization accuracy of an 
eight node ASLA formation is compared with an equal numbered circular formation at 
varying target bearing. From the simulation, we see that the bearing of the target has 
little to no effect on the circular estimate, whereas the linear formation performs better 


when the target bearing is approximately + 30 degrees from the array’s normal direction. 
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Even though conditionally optimal, its high accuracy is the key reason the linear 


formation was chosen for the proposed scheme. 
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Figure 7. Accuracy of position estimate for fixed linear formation (blue) and fixed 
circular formation (red). 


5. TDOA Outlier Measurement Diagnostics 


With the proposed localization technique implemented on an elevated, mobile 
sensor network, there will be increased measurement noise due to sensor position and 
other sensor network related errors. This increase in measurement noise creates a need 
for methods to mitigate the effects of statistical outliers or errors on parameter estimation. 
In this section, we cover some statistical tools used in our proposed localization scheme 


for robustness to such errors. 


To assess the influence of the i” measurement/observation in a least-squares 
solution, we can calculate the solution with and without it. Calculating the amount of 


divergence with and without the measurement is termed “single case diagnostics” [18]. 
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To perform this diagnostic, it is important to consider which measure of divergence is 
appropriate for the solution. A standard approach is to use the squared-Mahalanobis 


distance given by [18], [19], 
Dy; = (v, -E{v,})C," (v, -E{v,})" > (27) 


where the v,,is the x2 matrix containing M location estimates, v, is the 1x2 vector 
containing the i" location estimate, and C,is the covariance matrix of v,. The purpose 


of the Mahalanobis distance is to discern how much influence the i” measurement has on 


the solution given the covariance of each of the measurements. 


To illustrate this concept, we consider a set of 200 weighted least-squares location 
estimates shown in Figure 8. In this set, there is one outlier shown in green, which we 
wish to reject. Using the Euclidean distance from the mean as a measure of divergence, 
we see that the outlier has the same divergence as the standard estimate shown in red. In 
contrast, using the Mahalanobis distance as the measure of divergence, we find that the 
outlier has a divergence of 34.5 where as the standard estimate has a divergence of 0.9; 


hence, the outlier can be easily identified and discarded. 
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Figure 8. Scatter plot of 200 hyperbolic location estimates. 
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Le BEAMFORMING 


Beamforming is an efficient signal processing technique that enables space- 
division multiple access [6]. The benefits and various applications of traditional 
beamforming have been well documented [60], [61]. The extension of beamforming to 
wireless sensor networks has given rise to the concept of collaborative beamforming [2], 
[62] in which a network of synchronized and distributed sensors is used to form a 


distributed beamforming array. 


Beamforming begins with a plane wave signal impinging on the array from a 
direction-of-arrival 0, With the signal arriving at each node with a corresponding phase 
delay, beamforming uses phase shifters to nullify this delay and constructively sums each 
node’s signal. The concept of a plane wave falling across a uniform linear array of M 
nodes where the phase delay is a function of the distance between each node d is 


illustrated in Figure 9. 


Array normal 





Array node i= 1 2 BSB ——————aaeve M 


Figure 9. Plane wave impinging on a uniform linear array. 
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1. Phase Shift Beamforming 


With the phase of the incoming signal at node i=1 set to zero, the combined 
signal of a uniform linear array with M nodes, commonly termed the array factor can be 


expressed as [63] 


M 
4,66 )=y wer, (28) 


i=l 
where 


w, = Ae ere (29) 


is the complex weight with magnitude 4; and phase —fdsinO@,, A is the signal 
wavelength and f= (27) / 2. From (28), it is apparent that the array yields a maximum 
response in the direction of the beam steering angle 0. By controlling the response of 


the array as a function of @,,, the spatial diversity required for space-division multiple 


access is achieved. Along with this diversity, beamforming also provides a signal gain 


equal to the number of receiving elements/nodes. 


There are many parameters that dictate the resulting array factor of a 
beamforming array. Of key importance to this research is the number elements/nodes 
and their placement relative to each other. The effect of increasing the number of 
elements/nodes in a uniform linear array’s array factor is illustrated in Figure 10. The 
first effect is the narrowing of the main beam’s beamwidth. The main beam is defined as 
the lobe that contains the direction of the array factor’s maximum response. The 
beamwidth, also called the half-power beam width, is typically defined as the angle range 
at which the main beam falls above 3 dB of the maximum response. The relationship 


between the beamwidth and the number of nodes in a linear array is given by [63] 
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Bix 0.882 . 
dM 





(30) 


The second effect is the location and the number of sidelobes present in the array factor. 
Sidelobes are defined as any lobe other than the main beam lobe. For uniform linear 
arrays, the magnitude of the nearest sidelobe to the main lobe is on average 13 dB lower 


[63]. 


Array Geometry 





90 


Normalized Power (dB) 


Figure 10. Array factor of a 15 isotropic node linear array with d=1/2. 


2. Adaptive Beamforming 


Adaptive beamforming is the process of adapting an antenna’s array factor in a 
manner that both increase the gain of a target signal while decreasing the gain of all 
interfering signals. This is accomplished by placing a complex weight on each antenna 
signal. These weights modify the phase and allow for the coherent amplification 


(reduction) of a target signal (undesired signals). In this section, we will illustrate the 
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basic phase shift beamformer [64], [65] and the adaptive minimum variance distortionless 


response (MVDR) beamformer [66], [67]. 

The MVDR beamformer combats the effects of interfering signals by placing 
nulls in their direction. This is achieved by selecting the array weights that minimize the 
output noise variance [63]. The array output is expressed as 


y(t)=w"s+win, (31) 


where w is an Mx1 vector of the complex weights to be determined, u is an Mx1 


vector of the sum of all interfering signal vectors 7, =[l,e”“""’,...,e” eee 
represented as 
N, 
w= Di (t)M» (32) 
and s isan Mx\1 the signal vector represented by 
s=s(t)v, (33) 
with s(t) being the desired signal and the steering vector v =[l,e“""",...,e’” eae i 


By minimizing the noise contributions of the interfering signals in the output, we can 


determine the optimal weights w,. The output variance can be shown as [63] 
2_.H H 
o,=w Rwtw' Rw. (34) 


Assuming that R, = E {ss \ and R, =E { uu" \ are known, we see that the minimization 


of the output variance is reduced to the minimization of w”R,w, which results in an 


expression for the optimal weights [63] 
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Wi= az, Y. (35) 


To illustrate both the phase shift and MVDR beamforming concepts, we consider 
a scenario where a target signal with two interfering signals is received by the 15 node 
uniform linear array shown in Figure 10. The target signal-of-interest is shown in 
Figure 11, with a direction-of-arrival (DOA) of 45 deg and a signal amplitude of one. 
The two interfering signals are arriving at a DOA of 30 and 50 deg, each with an 


amplitude of 10. 
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Figure 11. Transmitted signal-of-interest. 


The phase shift beamformer’s response to the interfering signals and target signal 
is shown in Figure 12. We can see the two interfering signals have fully corrupted the 
target signal. In contrast to the phase shift beamformer, the adaptive MVDR beamformer 
response is shown in Figure 13. From this, we can clearly see that the MVDR 


beamformer is able to nullify interference and successfully recover the target signal. 


28 





a =e = = Po Po 
ie} oS oO foo) oO ho 
T T T T T T 

1 1 


Magnitude (V) 
=) 


0 0.05 0.1 0.15 0.2 0.25 0.3 
Time (seconds) 











Figure 12. Response of 15 node phase shift beamformer with two 10 dB interfering 
signals at DOA of 30 and 50 degrees. 
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Figure 13. Output signal of the MVDR beamformer in the presence of a target 
signal at 0. =45deg and two interfering signals at 9, =30 and 50 deg. 


A comparison between the array factor of the MVDR and that of the phase shift 
beamformer is provided in Figure 14. From the illustration, we can see the MVDR’s 


array factor contains two strong nulls in the directions of the interference, whereas the 
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phase shift beamformer has a response of approximately —13 dB and —3 dB in the same 


directions. 
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Figure 14. Comparison between the MVDR and phase shift beamformers in the 
presence of a target signal at 0 =45deg and two interfering signals at 


@. = 30 and 50 deg. 


3s Error Effects in Beamforming 


Two key errors to account for in beamforming are pointing errors and unintended 
sidelobe response. Pointing errors are defined as the difference between the array’s 
actual steering direction and that of the maximum response. The standard deviation of 


the pointing error is given by [4] 


(36) 


where o, is the standard deviation of the phase errors. With the beamwidth for a 


uniform linear array as defined in (30), the pointing errors can be expressed as 


30 


0.8840, 





oO ~———_, a 
PEULA dM JM (37) 
The expression for the expected side lobe level is given by [63] 
a(1-A*)+ox, 
A, = ( ; (38) 


Ma’ A’ 


2 : : a, x ‘i ‘ 
where a is the error free main beam gain, oj, is the variance of amplitude error and 


2 
Ong 


od z 
Ars le {e*} . If the phase errors are zero mean and Gaussian, then A” =e 


4. Sidelobe Control via Array Tapering 


For a uniform linear array, the sidelobe response can be controlled by adjusting 
the gain of each node’s response [68]. This process is termed tapering and is akin to 
windowing in digital filtering [69], [70]. Tapering an array provides a tradeoff between 
the beamwidth and the sidelobe levels. With this in mind, we will briefly examine the 
characteristics of four tapering methods: Uniform [71], Dolph-Chebyshev [72], Taylor- 
Kaiser [73], and Binomial [74]. 


a. Tapering Methods for Uniform Linear Arrays 


The uniform tapered array is actually a non-tapered array as all nodes are 
uniformly weighted by one, which results in the standard phase shift beamformer with an 


array factor of 


M 
An (0) A yore (39) 


m=1 


This tapering method makes no consideration to sidelobe response and has the highest 


sidelobe level when compared to other methods. 
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The binomial tapered array is unique because it has zero sidelobes, but this is at 
the cost of a wider beamwidth. The weights of an M-node binomial uniform linear array 


are calculated from the binomial coefficients as follows [70] 





m=0,1,...,M—I1. (40) 


The corresponding array factor is given by [70] 


drn(0)=[206{)] (41) 


From this expression, we can see that the array factor decreases to zero as @ approaches 
ae 70%. 

The Dolph-Chebyshev [72] taper method is best applied when the sidelobe 
response at all angles must be kept below a specific value Rp. The taper weights for the 


Dolph-Chebyshev uniform linear array are calculated using the Dolph-Chebyshev 
window transform [70], [75] 


{soe eof 








Were (@; ) = . i=O01,..,M-L., (42) 
cre (e) cosh (Mu cosh! (R, )) 
where the scaling factor R, is defined as [75] 
cosh! (10% ”) 
R, =cosh ; (43) 
M 


The array weights w,,.(m) are then computed as the inverse discrete Fourier transform of 


Wax (@,). The corresponding array factor for an array with an even number of nodes is 


expressed as [70] 
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App (9) = cos 2 cos’ [x oos(2))] (44) 


The Taylor-Kaiser tapering method is similar to the Dolph-Chebyshev method but 
has an exponentially decreasing sidelobe response. The array weights can be computed 


as [75] 


Wig (1m) =I (Atty Vim? M2), m=1,...,M (45) 


where J, is the zero”-order modified Bessel function of the first kind and @,, is the 
0 TK 


parameter used to control the sidelobe level response. The Taylor-Kaiser array factor can 


be expressed as [70] 


M 
1h (207 ) fa2,2° —(M6/2) 


sinh{ an a -(M6/2) 





(46) 





Arr (0) = 


b. Comparison of the Taper Methods 


A comparison of the four taper methods is illustrated in Figure 15. From the 
results, we see that the beamwidth is inversely proportional to the sidelobe response. We 
can see that the uniform tapering achieves the narrowest beamwidth with the highest 
sidelobe response, whereas the binomial tapering has the widest beamwidth with zero 
sidelobe response. The Dolph-Chebyshev has the unique ability to maintain a given 
maximum sidelobe response, while the Taylor-Kaiser taper show a better sidelobe 


response at the cost of a slightly wider beamwidth. 
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c) Taylor-Kaiser taper d) Binomial taper 


Figure 15. Beam width and sidelobe level tradeoff for four different tapering 


methods using a 15 node ULA. The beamwidth 3 dB |ine is shown in red, 
from [70]. 


=P Grating Lobe Control via Virtual Filling 


In collaborative beamforming where the nodes are separated by a distance greater 
than 4/2, grating lobes in the array factor will be a significant issue [71]. Grating lobes 
are lobes that have the same magnitude as the main lobe. To minimize the effects of 
grating lobes, the array can be virtually filled [76]. The virtual filling technique 


minimizes the grating lobes by filling in any internode distance d, greater than 2/2 
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with virtual nodes. These virtual nodes consist of interpolated complex signal data 
derived from the real nodes. This filling is done in a manner that turns the array into a 


virtually filled uniform linear array, which is inherently free of grating lobes [76]. 


A single snapshot of the complex signal received at the n" node in a collaborative 


beamforming network can be expressed as 
Nr 
A, ae exp(j¢,; )exp(-jBx, sin ,,), n=l...,.M, (47) 
i=l 


where JN, is the number of signals impinging on the array, x, is the n" node’s location 


n 


h 


on the x-axis, V, signal’s magnitude, phase, and bearing, 


ti? 


t 
g,, and @, are the i 


respectively. This expression can be written in matrix form as [76] 
(48) 


where A, isan Mx] vector containing the complex signal snapshot from each node, E, 


isan M x N, matrix given by 


exp(—j Bx, sin 6, ) ve exp(— jx, sin 8, y. ) 
E= ay ‘ (49) 
exp(-/Px,sin8.,) - exp(—jBx, sind, x, ) 
and the NV, xl vector 
V,, exp (4, ) 
Vie (50) 
Viv, &XP (id, ) 


Assuming that each signal’s bearing is known, an estimate of their magnitude and phase 


can be obtained using the least-squares approach expressed as [77] 
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V, =(ETE,) ETA,. (51) 


These signal estimates are then used to represent each virtual nodes signal output. The 


output of each virtual node is expressed as 
Np A nN 
A, =>, ,exp( Jd, Jexp(-i8x, sin9,,), v= 1.4), (52) 
i=l 


where x. is the position of the v" virtual node along the x-axis and M. is the number of 
\ p g v 


virtual nodes equal to d,/(A/2). 


In this chapter, we covered hyperbolic localization and beamforming within the 
context of wireless sensor networks. With regard to localization, we presented two 
localization techniques and their application to an optimal ASLA sensor formation. We 
also covered the use of outlier detection techniques to mitigate measurement noise. The 
discussion regarding beamforming focused on techniques to increase array robustness 
against various steering errors and sidelobe interference. Together, these topics provide 


the necessary details to present our proposed signal collection scheme in the next chapter. 
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ll. SOLUTION APPROACH 


This chapter begins with an explanation of the operating scenario and 
assumptions used in this dissertation. We provide an overview of the proposed scheme 
for radio frequency source localization and signal collection. We then define the figures 


of merit used to evaluate the scheme’s performance. 


A. PROPOSED SCHEME 


The intended operational scenario for the proposed scheme is shown in Figure 16. 


With a distant, stationary, sporadic or repeating source emitter located at (x,,y,), there is 


a network of sensors deployed with the goal of locating the emitter and collecting its 
signals. In this scenario and in this dissertation research in general, we assume the 
following. Sensor node deployment is restricted to a circular area, each node has an 
isotropic antenna with matched polarization, and all nodes are synchronized to a common 


clock. 


The network of M sensor nodes is initially deployed into an ASLA formation 
(introduced in Chapter II) where they form an ad hoc wireless sensor network. Each 
node consists of a multirotor UAV, used to realize the concept of an elevated, mobile 
sensor network. The ASLA formation shown in Figure 6 places the M sensor nodes into 
three node stacks, each containing an equal number of nodes. These three stacks are 
deployed along the x-axis, one at the origin and the other two flanking the center on both 


sides at a distance of 7,. For source localization, this formation is used to minimize the 


estimate error variance [17], [34]. For signal collection, the node stacks allow for robust 


signal estimation. 
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Figure 16. Operational concept of an elevated, mobile sensor network deployed in 
an ASLA formation. 


A schematic diagram of the proposed source localization and signal collection 
scheme is shown in Figure 17. The scheme begins in the localization phase once a signal 
has been detected. This generates a set of TDOA measurements that is fed to the location 
estimator. Using this initial location estimate, we reorient the sensor network to be 
perpendicular to the target emitter. The objective of this network orientation is to 
maximize the accuracy of the location estimate, which we call the refined estimate. With 
the refined estimate obtaining an improved location estimate, this information is used in 
the signal collection phase. In this phase, each node in the network samples the signal 
and transmits these samples to the network’s sink node. Located in the center of the 
network, the sink node uses a combination of beamform processing and signal estimation 
to combine and amplify the signal samples coherently. This makes the signal collection 


phase essentially a collaborative beamforming effort. 
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Figure 17. Process flow for propose signal collection scheme. 


1. Source Localization 


The objective of this phase is to obtain a precise location estimate of the target 
signal emitter. To accomplish this, we employ the use of an unbiased and efficient 
maximum-likelihood location estimator [52]. Although an effective location estimator, it 
requires an initial position estimate, in order to avoid local minima solutions [59]. To 
satisfy this requirement, we precede the maximum-likelihood estimator with another 
estimator that does not require initialization. For this task, we propose a weighted least- 


squares estimator. 


Furthermore, after obtaining the initial estimate we take advantage of the 
network’s mobility and re-orient the network nodes to create an optimal sensor to target 
geometry [17]. This new orientation minimizes the geometric dilution of precision [52], 
thus minimizing the maximum-likelihood estimator’s error variance. Because of this 
further refinement of the location estimate, we term the maximum-likelihood estimator 
the refining estimator. This two-stage localization technique is illustrated in Figure 18. 


Starting in the initial configuration shown in Figure 18 (a), we obtain a location estimate 
(%,,9,) which is used to calculate the emitter bearing estimate 6. Using 6,, we reorient 


the network nodes perpendicular to the target signal emitter, resulting in the reoriented 


network configuration as seen in Figure 18 (b). 
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Figure 18. Network reorientation operation for the two-stage localization technique. 


In conjunction with both estimators, we propose to implement a robust 
measurement outlier rejection process to detect and reject erroneous TDOA 
measurements. This process is used to increase the localization phase’s robustness to 


measurement error and sensor position errors. 


2. Sensor Node Position Error 


In this research, the position information of each node is assumed known. 
Fluctuations in each node’s position due to their station keeping operations are also 
assumed [47], [48], [78]. The effects of these position errors are the main concern in this 
dissertation. With each node in the network envisioned as a multirotor UAV, the actual 
position errors are mainly governed by wind, GPS accuracy, and flight controller scheme 
[46]-[48]. In the absence of such knowledge we resort to two basic position error 
distributions, uniform and Gaussian. More specifically, we model these small 


fluctuations in the x and y positions as two independent and identically distributed 


uniform or Gaussian random variables 6, and 6,,, respectively. 


For the uniformly distributed position error case, the probability density functions 


are given by 
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tx (8-55 2055020. 
(53) 
fe. esa -6, <6, <6, 


where 6, can be interpreted as the maximum position error in the x or y direction [48]. 
For the case of Gaussian position error, both random variables are zero mean with 


a variance of om The probability density functions are expressed as 








oO 
FS, (6,)= exp ols -0 <6, <0 
: ( ) 2n0;, 20, 
. (54) 
1 
Je, |= exp 7_|, —0 <6, <0 
‘ ( ) 2x07, 26, . 


where ©”, is the position error variance in the x or y direction [48]. 


3. Collaborative Signal Collection 


After the source localization, the objective of the signal collection is to collect and 
amplify the target signal. The main concern is that, while the mobility of each node is 
allowed for optimal formations in the localization phase, it will cause array phase errors. 
To achieve the objective, despite the presence of such errors, we propose a combination 


of collaborative beamforming and signal estimation. 


In this approach, we exploit the sensor grouping of the ASLA formation. The 
ASLA formation contains three groups of sensors placed symmetrically along the x-axis 
(see Figure 6). With each node in a group providing a noisy sample of the same signal, a 
sample mean can be calculated. Using each group’s sample mean in a collaborative 


phase shift beamformer, we reconstruct and amplify the signal. 
Furthermore, since the ASLA formation’s array factor contains grating lobes due 
to its large inter-stack distancer,, it is extra sensitive to interfering signals. To increase 


the robustness of signal collection against such signals, we employ the use of virtual 
4] 


filling [76]. This array processing technique is used to manage the array’s inherent side 


lobes and grating lobes in order to minimize the effects of interfering signals. 


It is important to note that while our scheme is not a traditional adaptive 
beamformer, it accomplishes the same task. Adaptive beamforming is the process of 
manipulating an array’s weights in order to amplify a target signal while suppressing 
undesirable signals. Although our approach uses fixed uniform weights, it amplifies a 
target signal using array rotation and suppresses undesirable signal through the use of 


virtual filling and tapering techniques. 


B. PERFORMANCE METRICS 


With this dissertation heavily focused on parameter estimation, it is critical to 
define performance metrics as a means to support and analyze the proposed theory and 


scheme. The performance metrics used throughout this dissertation is defined as follows. 


For localization schemes, a widely used performance metric is the Cramer—Rao 
lower bound (CRLB) [58]. The CRLB is the theoretical minimum solution variance an 
estimator can achieve. Any unbiased estimator that can achieve the CRLB is said to be 
efficient [59]. For our proposed network, the CRLB for location (minimum location 


uncertainty) and bearing (see Chapter II) are expressed as 





-1 
Mr* : 
tua"| is sin’ | (55) 
R t 
and 
2Mr° wo ip 
cere = Boe sin’ 8, |; (56) 
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where M is the number of sensor, 7, is the distance between the outer and center node 

stacks, oj, is the ranging measurement error variance, and 6, is the target source bearing. 
For estimation processes, the root mean-square error is a standard performance 

metric. For source localization, we focus on two root mean-square error values, the first 


being root mean-square error of location estimates defined as 
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where @, is a 2x1 vector containing the true source’s location, @, is the estimate @,, and 
n,is the number of trials. The second value is the root mean-square error of bearing 


estimates defined as 


(58) 





where @, is the true source bearing and 6 is its estimate. 


For signal estimation, we focus on the root mean-square error of magnitude and 


phase. The root mean-square error of magnitude estimates is defined as 


(59) 





where V, is the true signal magnitude and V, is its estimate. The root mean-square error of 


phase estimates is defined as 


(60) 





where @ is the true signal magnitude and d, is its estimate. 


To quantify the error in a normalized array factor with and without position 


errors, we define the normalized array factor error as 


v.=["(4(8)-4,(8)) 49, (61) 


where @ is the signal bearing, A, is the normalized error free array factor, and A, is the 
normalized array factor with position errors. 


In this chapter, we provided an overview of the proposed scheme for source 
localization and signal collection. We described the use of two sequential location 
estimates that support a follow on signal collection. We then described a set of 
performance metrics used to analyze the scheme’s performance. In the next chapter, we 
present a two-stage localization technique capable of approaching the CRLB and propose 


a measurement outlier rejection process to increase localization robustness. 


4A 


IV. ROBUST TWO-STAGE SOURCE LOCALIZATION FROM AN 
AIRBORNE SYMMETRIC LINE ARRAY NETWORK 


In this chapter, we introduce a two-stage source localization technique. We 
examine its performance in the presence of measurement and sensor position errors. To 
combat these effects, we propose a measurement outlier rejection technique. This 
technique is used to identify and reject specious TDOA measurements in an effort to 
increase the localization robustness. We then use simulations to support the theoretical 


development and analyze the performance. 


A. TWO-STAGE HYPERBOLIC LOCALIZATION 


The objective is to estimate the location of a single source emitter located at 


(x,,y,)- The sensor network is deployed in the initial ASLA configuration, as seen in 


Figure 18 (a). In this configuration, the network detects an incoming emission and 
generates a set of TDOA measurements. Using these measurements in the weighted 
least-squares estimator introduced in Chapter II, we obtain a location estimate. This 
initial localization is a closed-form estimator that does not require any initialization [5], 


[79]. 


1. Hyperbolic Localization via Weighted Least-Squares 


The estimator described in this section is derived from the linear array estimators 
found in [5] and [79], and a generic linear array formation is adapted for use with an 


ASLA formation. 


Since weighted least-squares estimation is primarily for linear models, its 
applicability to hyperbolic estimation is not immediately obvious. To apply the weighted 
least-squares approach to localization, we must rewrite (6) to obtain a set of linear 
equations. To do this, we define the noisy range difference-of-arrival measurement of a 


sensor pair as [5], [79] 


di. =a, +E 5, i,k 7 eee if ; (62) 
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0 ; : ‘ : 
where d;, is the true range difference between sensors i and k and é€,, is the zero mean 


Gaussian range difference measurement error. Then the range from the i" sensor to the 


source emitter is determined as [5] 





nav(u-x) t0r-y) - (63) 
By squaring both sides of (63), we obtain 


rt =K,-2x,x,-2y,y,+x7 +y, i=1,...,M, (64) 


where K,=x>+y? and r,,=cd,,=",-7,. In this form, we can formulate a linear 
weighted least-squares problem in which the sensors are placed in a line. By substituting 


2 2. 
1? =(r,, +7) in (64), we get 


my +2rn nth =K,-2x,x,-2y,y,+K,. (65) 


By co-locating the origin with the first sensor and making it the reference sensor, we get 
r,=r,,and K, =K,, then (65) can be expressed as 
i +24 =—2X,x, -2y,,y, + K, —K,. (66) 


1, 


Finally, with all the sensors in a line, we replace —2x,,x,—2y,,y, by —2x,, (x, +a, y,) 


where @, is some constant. With this substitution, (66) becomes 


= —K. + K, — 2%, 4 (% F: ay, ) = 2h 7 (67) 


i; i 
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With (67) now linear in 7 and (x,+a,y,), a weighted least-squares solution can be 


obtained as [5], [79] 


(x +a y ) -1 
A = t wt) | (ATp-i P| 
Ors -| : =(G'C,'G,) GiC,'h, (68) 
t, 
where 
7 2 2 
Xo, Ky, —X; 
2 2 
Xx iz Ky SX, 
| at 31 _| 817% 
G, aed : p) h =< p) 
2 2 
Xu lui 1 *u 





and Cc, is the inverse covariance matrix of the TDOA measurement vector from (62). 


When all the sensors are located on the x-axis, y,,=Ofor i=2,...,M, resulting in 


a, =9. The solution is then expressed as 


Ons = | 


Using (69), we calculate the y-coordinate estimate using the expression 


> 


> 


( =(G'C3G,) GICSh. (69) 


~ 


jh, =P? -x2. (70) 


Using the initial estimate @,,,, we see that the network reorients itself into a 
desired ASLA formation, as seen in Figure 18 (b). In this configuration, the effect of 
GDOP is minimized, thus minimizing the estimate’s CRLB [34]. To take full advantage 
of the new formation and its minimized CRLB, we propose the use of a maximum- 
likelihood estimator introduced in Chapter II. This estimator is unbiased and 


asymptotically efficient, 1.e., it achieves the CRLB [59]. 
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2, Hyperbolic Localization via Maximum-Likelihood Estimation 


The estimator described in this section is derived from the generic TDOA 
estimator found in [52] but here is adapted for use with an ASLA formation in polar 


coordinates. 


Introduced in Chapter II, the maximum-likelihood estimator has been shown to be 
unbiased and asymptotically efficient [59]. Its extension to hyperbolic localization is 
done in a similar fashion to that of the weighted least-squares estimator. We rewrite (13) 


to account for time-of-arrival measurements as 


t=t)1,+D/cteE (71) 


where D is an M x1 vector containing the range from the emitter to each node, f, is the 
time of signal emission, J,, is an Mx1 vector of ones, and ¢ is an Mx\|vector 


containing the TOA measurement errors. To convert the TOAs to TDOAs, we subtract 
the 7" TOA measurement from the first TOA measurement to eliminate t,. This results in 


t,-t, =t-4) +(D,-D,)/c+4,-6,. (72) 


To put it in matrix form, we multiply (71) by the (/—1)xM matrix [52] 


1 -1l 0 0 
rome 10-1. : (73) 
we lay es ae ae 0 
1 0 0-1 
to get 
Vip =GyD/ct+Gye. (74) 
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where y,, =G,,t is an (M-1)x1 vector containing all the resulting time-difference-of- 


arrival measurements. 


From (18), we set @. =@,,., and the maximum-likelihood estimator’s solution in 


polar coordinates can be expressed as [52] 


Ou -| 


where @,,; is the weighted least-squares estimate in polar coordinates, D, is the 


é 


j Bn te(M GCG HY HCE (x9 (G,D,)/*). (75) 


(M -1)x1 vector containing the distances between the reference point and the sensor 


nodes and H is the (M -1)x2 matrix expressed as 








| r,,-r.cos(6,-6,) — r,r,sin(@,, -4,) 
Dix Dy. 
H= ‘ : ; m=2,...,M, (76) 
ry —1,cos(O,—-9,) nyr,sin(O, —9, ) 
Dou Dou 





with 7, and @, the m™ 


node’s range and bearing, respectively, 7. and @. are the 
reference point’s range and bearing, respectively, D,,, the m'" element of D,, and the 


product G,,H the hyperbolic version of (17). 


B. LOCALIZATION PERFORMANCE 


A series of Monte Carlo simulations to support proposed localization technique is 
presented in this section. The results of the simulation are based on 10,000 trials each. 
For each simulation, the ASLA configuration is arranged in accordance with Figure 18, 
with 7, = 200 m and NV, =5. Note, since TDOA measurements are typically on the order 
of nanoseconds, for convenience, we use range difference-of-arrival (RDOA), i.e., the 


TDOA measurement weighted by the signal propagation speed c, which is typically on 
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the order of meters. This substitution is used only as a means to simplify the simulation 


results and does not affect the theory in any way. 
In Figure 19, the root mean-square error of location estimate ¢,,. for both the 


initial and refined localization are compared against their corresponding CRLB as a 


function of RDOA noise variance 0%, =c’o;,. In this simulation, N, =5, the emitter is 
located at w, =[2500 m, 45 deg], r,=200m, and each data point is the result of 
10,000 trials. From the results, we see that the ¢,,. of both the initial and refined 


estimates increase rapidly with noise. As expected, the ¢,,. of the refined estimate 


outperforms the initial estimate at all values of noise simulated. Also, we can see that the 
refined estimator using a maximum-likelihood estimate approaches the CRLB up to 
approximately 0 dB. It is important to note that the maximum-likelihood estimate is 
asymptotically efficient [15], but here, using the ASLA formation, it approaches the 


CRLB in one iteration. 
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Figure 19. Root mean-square error ¢,,. of the initial (black) and refined (red) 
localizations versus RDOA noise. Result based on 10,000 trials, NV, =5, 
r, =200 m, and @, =[2500 m, 45 deg]’. 
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A plot of the angular root mean-square error ¢, for both the initial and refined 
localization as a function of Os is provided in Figure 20. In this simulation, V, =5, the 
emitter is located at w, =[2500 m, 45 deg]", 7, = 200 m, and each data point is the result 


of 10,000 trials. The y-axis is €, in radians, and the x-axis is RDOA noise in dB. From 


the results, we observe that the same trends as in Figure 19. From (56), we know that the 


ASLA formation yields an accurate bearing estimate. This is validated by the simulation 
results with the refined estimator’s €, =—0.25x10~ rad at -10 dB RDOA noise. The 
initial bearing estimate is just shy of CRLB, while the refined maximum-likelihood 
estimate approaches it at all simulated values. As with the ¢,,. estimate, the maximum- 


likelihood estimator also approaches the CRLB in one iteration. 
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Figure 20. RMSE €¢, of the initial (black) and refined (red) localizations versus 
RDOA noise. Result based on 10,000 trials, VN, =5, 7, =200 m, and 
cw, =[2500 m, 45 deg]'. 
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We plot both the initial and refined localization’s ¢,,. values in Figure 21 as a 
function of sensors per stackN,. Each estimator’s CRLB is also included for 
comparison. In this simulation, o;%,=0dB, the emitter is located at 
o, =[2000 m, 30 deg]', r,=200 m, and each data point is the result of 10,000 trials. 
From the results, we see that the ¢,,. of the initial estimate approaches the CRLB at 
N, =2 but quickly diverges from it. In contrast, the ¢,,. of the refined estimate does not 


approach the CRLB until V, >5. It is interesting to note that an increasing N, improves 


Ss 


the performance of the refined estimate at a higher rate than the initial estimate. 
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Figure 21. Root mean-square error ¢,,. of the initial (black) and refined (red) 
localizations versus N,. Results based on 10,000 trials with o;,, =0 dB, 
N, =5, r,=200 m, and @, =[2000 m, 30 deg]’. 


A comparison of both the initial and refined localization’s ¢, values is shown in 


Figure 22 as a function of sensors per stack N,. In this simulation, o%, =0 dB, the 
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emitter is located at w, =[2000 m, 30 deg]’, r, = 200 m, and each data point is the result 
of 10,000 trials. From the results, we see that the ¢, of the initial estimate is always just 


shy of the CRLB for all values of N, simulated. Conversely, the ¢, of the refined 


Ss 


estimate approaches the CRLB for all values of NV, simulated. These results suggest that 


the refined maximum-likelihood estimator is a good candidate for bearing estimates 


regardless of the number of sensors employed. 
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Figure 22. Root mean-square error &, of the initial (black) and refined (red) 
localizations versus N,. Results based on 10,000 trials with o;,, =0 dB, 
N, =5, r,=200 m, and @, =[2000 m, 30 deg]". 


In Figure 23, both the initial and refined localization’s ¢,. are plotted as a 
function of the emitter’s bearing. In this simulation, Ors =0 dB, the emitter is located at 


o, =[2000 m, 6)" , ©, =200 m, and each data point is the result of 10,000 trials. From 


the results, it is clear that the initial and refined estimator’s performance is dependent on 


the target’s bearing (relative to the initial array configuration). The initial estimate 


a3 


follows the CRLB as seen in (55), while the refined estimate improves the ¢,,,. within the 
range of 40<@ <40deg. From (55), we see that ¢,. approaches infinity as 0 


approaches + 90 deg. This suggests that the initial estimate requires some knowledge of 
the target bearing a priori to insure a sufficiently accurate initialization for the refined 
estimate. In other words, if the emitter bearing is greater than + 40 deg, the initial 
estimate is insufficient in setting up the follow-on refined estimate. If this is not possible, 


the initial estimate can be obtained using estimators that are independent of @ such as 


the ones in [5] and [80]. Note that these estimators require a random positioning of the 


sensor nodes and are not applicable to an ASLA configuration. 
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Figure 23. Root mean-square error ¢,,. of the initial (0) and refined (*) 
localizations versus 0,. Results based on 10,000 trials with 0%, =0 dB, 
N,=5, 7, =200 m, andv, = 2000 m. 


1. Performance of the Localization Scheme in the Presence of Uniform 
Position Errors 


This research envisions the use of multirotor UAVs to realize an elevated, mobile 


WSN. As such, it is necessary to address the position errors induced by the UAV’s 
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unsteady station-keeping operations. This unsteadiness is mainly due to wind [46] and 
causes nodes to momentarily deviate from their given positions. In this section, we 
model these position errors as two uniformly distributed random variables 6, and 6,, 
each bounded by 6,. The true distribution of position errors for any given UAV is 


determined by many factors such as GPS accuracy, wind effect, and flight controller [38]. 
In the absence of such knowledge, we resort to a basic distribution in order to complete 
our analysis. The effects of these position errors are illustrated in Figure 24 through 


Figure 27. 


A comparison of the initial localization’s root mean squared location error ¢,,. 
and bearing error €, with and without position errors is shown in Figure 24 and 
Figure 25 as a function of the maximum position error é, (see (53)). In this simulation, 
Oo; =0 dB, the emitter is located at @, =[2000 m, 30 deg]", r, = 200 m, and each data 


point is the result of 10,000 trials. From the plots, we see that the root mean squared 


location error monotonically increases with 6,,, both in location and bearing estimates. 








—- No Position Errors 
—©- With Position Errors 























Figure 24. Initial localization’s root mean-square error ¢,, with no position errors 
(black) and with position errors (red) versus maximum uniform position 
error é, . Results based on 10,000 trials with Son =0 dB, N, =5S, 


r, = 200 m, anda, =[2000 m, 30 deg]’. 
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Figure 25. Initial localization’s root mean-square error &, with no position errors 
(black) and with position errors (red) versus maximum uniform position 
error 0,,. Results based on 10,000 trials with a =0 dB, NV, =5, 


r, = 200 m, anda, =[2000 m, 30 deg]". 


Similar to Figure 24 and Figure 25, the performance of the refined localization as 


a function of 6, is shown in Figure 26 and Figure 27. In this simulation, a. = 0 dB, the 
emitter is located at w, =[2000 m, 30 deg]', r, = 200 m, and each data point is the result 
of 10,000 trials. From the results, we see that the refined maximum-likelihood estimate 
is more resilient against position error with an almost linear response to 6,. Specifically, 
we see in Figure 26 that for é, =Im, ¢,,.. only increases by approximately 50 m. 


Overall, these results suggest that position errors play a significant role in the accuracy of 
localization using a mobile network. In the next section, we address these errors using a 


measurement outlier rejection process. 
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Figure 26. Refined localization’s root mean-square error ¢,,, with no position 


errors (black) and with position errors (red) versus maximum uniform 
position error 6,. Results based on 10,000 trials with Orn =0 dB, N, =5, 


r, =200 m, anda, =[2000 m, 30 deg]’. 
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Figure 27. Refined localization’s root mean-square error ¢, with no position errors 


(black) and with position errors (red) versus maximum uniform position 


error 0,. Results based on 10,000 trials with o,,=0dB, NV, =5. 


r, = 200 m, anda, =[2000 m, 30 deg]'. 
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C; ROBUST SIGNAL LOCALIZATION VIA MEASUREMENT OUTLIER 
REJECTION 


With the goal of improving the robustness of the source localization technique, 
we propose a measurement outlier rejection process to improve localization noise 
performance. This includes resilience to increased signal noise, sensor position errors, 
and measurement outliers. These improvements leverage the existing outlier detection 
techniques to identify and remove specious measurement data, i.e., measurements that 
seem benign but in actuality degrade the estimate’s accuracy. The basic premise is to 
determine the effect each measurement has on the end solution as a means to identify and 
remove measurement outliers. In this research, we define a measurement outlier as a 
measurement taken with a variance greater than a specified multiple of the expected 


variance. 


1. Robust Hyperbolic Localization via Measurement Outlier Rejection 


When trying to detect the presence of outliers in a data set containing location 
estimates, the simplest methodology is to calculate each sample’s Euclidean distance 
from the sample mean. The problem with this approach is that outliers heavily influence 
the sample mean, making it a poor discriminator of outliers [18]. Another approach is to 
use a statistical hypothesis test to determine the likelihood that the data set is derived 
from a single distribution [55]; however, this approach is more suited to problems with 


larger data sets (> 30 samples per stack) [55], which is not appropriate for our scheme. 


a. Measurement Outlier Rejection via the Mahalanobis Distance 


Since the localization technique is focused on a small sample size, the use of 
hypothesis testing to identify outliers is limited. An alternative approach is the 
Mahalanobis distance. This approach performs well with small sample sizes and is less 
affected by the presence of outliers [18]. In contrast to the Euclidean distance, the 
Mahalanobis distance calculates the distance from the mean but does so in a manner that 


accounts for the samples distribution, i.e., its covariance matrix. Given N,. location 


estimates, the corresponding Mahalanobis distance is defined as [18] 
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dy =(¥; -E{v,})C," (v, -E{v,}), (77) 


where the vis the N,.x2 matrix containing all the location estimates, v, is the 1x2 
vector containing the i" location estimate, and C,is the 2x2 covariance matrix of v,. 


In the proposed rejection process, we use the Mahalanobis distance of the weighted least- 
squares solution to determine the best and worst candidate measurements for use in the 
final location estimate. In this approach, we use single case diagnostics [18] to derive 
multiple weighted least-squares estimates from one set of TDOA measurements. For 


ASLA configured network, the number of single cases for a given signal measurement is 


Nye =M -1. 


A classic tool for the detection of statistical outliers in least-square estimation 
problems, single case diagnostics provide a methodology to single out measurements that 
lead to strong deviations in the end solution [18]. In this approach, a single measurement 
from a set is removed, and the corresponding weighted least-squares solution is 
calculated. This process is repeated until each measurement is treated in the same 
fashion. The Mahalanobis distance of each solution is then calculated using the mean 


and covariance of the entire N,. set of solutions. The solutions with high Mahalanobis 


distance values suggest that an influential measurement was removed. Such 
measurements are termed leverage points [18]. Conversely, if the Mahalanobis distance 


value was low, then the measurement in question is considered valid. 


b. Distribution of the Mahalanobis Distance 


With a measure of solution divergence now obtained, the next step is to be able to 
discern which values of this divergence represent outliers and which represent valid 
measurements. This process amounts to setting a Mahalanobis distance value outlier 
threshold. To calculate a suitable threshold value, we consider the distribution of the 
Mahalanobis distance values with and without outliers. Given a set of location estimates 
that are Gaussian with a known mean and covariance matrix, the distribution of the 


Mahalanobis distance values follow a Chi-squared distribution [81]. The Chi-squared 
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distribution is the resulting distribution for the sum of multiple squared Gaussian random 
variables. The Chi-squared distribution can be expressed by its degrees of freedom p, 
which is the number of squared Gaussian random variables present in the sum. The Chi- 


squared distribution of the Mahalanobis distance is expressed as [55] 


p/2-1 





1 
So, (du )= SHE (p72) P| dipt Nd) for 0<d,, <©, (78) 

where I’ is the gamma function defined as 

[(r)=[x"'ede — for r>0. (79) 

0 
The cumulative distribution function is defined as 
dy aA p-2)/2 4—x/2 
is (diz |?) = J 2E(p/2y (80) 


and the corresponding means and variance are expressed as [55] 


Ee dy} =D; 
us (81) 
oy, = 2p. 
Given that our scenario considers a small number of outliers, we hypothesize that the 


distribution of the Mahalanobis distance with and without outliers is significantly similar. 


To validate this hypothesis, a comparison of Mahalanobis distance values and 
their estimated distribution are shown in Figure 28. In this simulation, using (69), the 
histogram of 100 Mahalanobis distance values with and without outliers is shown. The 
Chi-squared probability density functions were plotted using the sample mean. In the 
outlier case, there are 15 outliers present with a variance five times larger than the 


expected variance. In both cases, the estimated degrees of freedom was p =1.98. These 
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results help support the claim that the estimated probability density function of the 


Mahalanobis distance approximates the Chi-squared with or without a small amount of 


outliers present. 
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Figure 28. Comparison of histogram and Chi-squared probability density function 


Cc. 


of Mahalanobis distances. 


Outlier Rejection Threshold 


Using the previous results to support our hypothesis, we next determine the 


outlier rejection threshold. A standard criterion for outlier rejection thresholds is the 


Tukey criteria for anomalies [82]. This criterion states that in a given range of data 


points, values at the extreme ends of the data set are considered outliers. To quantify the 


threshold, Tukey [82] proposed the following criterion. Taking the difference between 


the 25" and 75" percentile of a data set as the interquartile range [82] 


Qr = Q4 5% ~ Qys0 (82) 
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we consider any data point greater than Q.,.,+ (3/2) Q;p as an outlier [82]. It is important 


to note that this criterion makes no consideration on the data’s distribution. To 
incorporate knowledge of the data set distribution, we propose the use of Horn’s 


algorithm [83]. 


Horn’s algorithm begins with the estimation of a data set’s distribution. Next, the 
Tukey criterion threshold is determined using the estimated cumulative distribution 
function. Then, all data points greater than the calculated threshold are removed from the 
data set [83]. Overall, assuming the distribution of the Mahalanobis distance values is 


Chi-squared, the interquartile range is determined by [83] 


Qin =F '(0.75|p)-F'(0.25|p), (83) 


where the F'is the inverse Chi-squared cumulative distribution function. The outlier 


threshold using Horn’s algorithm can then be expressed as [82] [83] 


T, 


Oe = F'(0.75| p)+(3/2) Qe. (84) 
Using the outlier rejection process in conjunction with our proposed localization, 
we wish to minimize the effects of TDOA measurement noise, position errors, and 


measurement outliers. We term an estimator using this process a robust estimator. 


2. Performance of the Robust Localization in the Presence of 
Measurement Outliers 


To support the proposed outlier rejection process, we examine the results of four 
Monte Carlo simulations. Each simulation evaluates the robust localization estimator’s 
performance in the presence of random measurement outliers. Each simulation uses a 15 


node ASLA_ configured network in which the emitter is located at 


w, =[2000 m, 30 deg]’. 
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A comparison of the robust and regular weighted least-squares estimator’s root 


mean-square error €,, as a function of 0%, =c’o;, is shown in Figure 29. In this 
simulation, there are no measurement outliers present, N, =5, the emitter is located at 
o, =[2000 m, 30 deg], r,=200 m, and each data point is the result of 10,000 trials. 


The results show that the robust estimator’s ¢,,.is only slightly higher than the regular 


estimator at low noise levels. In contrast, the robust estimator outperforms the regular 
estimator at high noise levels. Overall, this result validates the effectiveness of the 


rejection process against high measurement noise. 
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Figure 29. Root mean-square error ¢,,. of the weighted least-squares (red) and 
robust weighted least-squares (blue) versus RDOA noise variance. No 
outliers. Results based on 10,000 trials with NV, =5, r, =200 m, and 


, =[2000 m, 30 deg]’. 


In Figure 30, the robust weighted least-squares estimator’s ¢,. is compared 
against its corresponding CRLB at increasing RDOA noise values. In this simulation, 
there are two random measurement outliers with a RDOA noise variance that is five 
times greater than the remaining set of measurements. Also, N,=5, the emitter is 
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located at w, =[2000 m, 30 deg]’, r, = 200 m, and each data point is the result of 10,000 


trials. The y-axis is ¢,,. in meters, and the x-axis is RDOA noise variance in dB. The 


results show that the effects of the outliers are almost negligible at the lower levels of 
RDOA noise (< —5 dB). The robust weighted least-squares estimator outperforms the 
weighted least-squares estimate at all simulated noise levels, with significant 
improvements occurring at RDOA noise level above —2 dB. Of interest is the robust 
estimator’s performance above 6 dB, where it again outperforms the no outlier case. This 
result validates the effectiveness of the proposed outlier rejection process for a wide 


range of measurement noise, with or without outliers present. 
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Figure 30. Root mean-square error ¢,,. of the weighted least-squares estimator with 


outliers (red), weighted least-squares estimator without outliers (black), and 
robust weighted least-squares estimator with outliers (blue) versus RDOA 
noise variance. Results based on 10,000 trials with NV, =5, r,=200 m, 


and , =[2000 m, 30 deg]'. In the outlier cases there were two random 


outliers present at a RDOA noise variance five times o;,. 
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The results of two Monte Carlo simulations illustrating the robust weighted least- 
squares estimator’s performance as a function of the number of random outliers present 


and their variance are shown in Figure 31 and Figure 32. 


A plot of the robust and regular weighted least-squares estimator ¢,,. is shown in 
Figure 31 as a function of the number of random outliers present. In this simulation, 
N,=5, the emitter is located at @, =[2000 m, 30 deg]’, 7, =200 m, o;,, =0 dB, the 
outlier variance is five times larger than the expected variance, and each data point is the 
result of 10,000 trials. From the results, we see that the increase in ¢,,. of both the 


regular and robust estimate is almost linear, with the robust estimator having a slower 
rate of increase. This result suggests that the performance increase of the robust 
estimator is a function of the number of outliers present. As the number of outliers 


increases, the more of them are rejected. 
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Figure 31. Root mean-square error ¢,,. of the weighted least-squares estimator with 


outliers (red), weighted least-squares estimator without outliers (black), and 
robust weighted least-squares estimator with outliers (blue) versus number 
of outliers. Results based on 10,000 trials with N, =5, r, =200 m, and 


@, =[2000 m, 30 deg]". In the outlier cases their RDOA noise variance is 


: 2 
five times o;,. 
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A comparison of the robust and regular weighted least-squares estimator’s ¢,,.. 


values as a function of outlier variance gain is shown in Figure 32. In this simulation, the 


parameters are similar to the ones used to create Figure 31. From the results, we see that 
the increase in ¢,,. approaches a maximum at variance multiplier greater than ten. Since 
each outlier is derived with a larger variance, they are easier to identify and be removed, 
thus we see a plateauing effect on the ¢,.. This also suggests that above a certain 
variance, value the ¢,,. of the robust estimator has a stronger dependence on the number 


of outliers rather than on their variance. 








—— WLS with No Outliers 
—=>— WLS with Outliers 
—>— Robust WLS with Outliers 











Sloe (M) 
© 
oO 








1 Ll L Ll i 1 Ll 1 i 1 
" 2 4 6 8 10 12 14 16 18 


Oultier Measurements RDOA noise variance (055) 





Figure 32. Root mean-square error ¢,,. of the weighted least-squares estimator with 


outliers (red), weighted least-squares estimator without outliers (black), and 
robust weighted least-squares estimator with outliers (blue) versus RDOA 


noise gain. Results based on 10,000 trials with NV, =5, 7, =200 m, and 


o, =[2000 m, 30 deg]'. In the outlier cases there were 2 random outliers 
present. 


The performance of the robust weighted least-squares estimator compared to the 
regular weighted least-squares estimator as a function of the emitter’s bearing is shown in 


Figure 33. In this simulation, there are two outlier measurements present with a variance 
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equal to five times o;,,, the emitter is located at w, =[2000 m, O]' where @, ranges 
from —60 <6, <60, 7, =200 m, and each data point is the result of 10,000 trials. From 
the results, we see the robust weighted least-squares estimator performs better than the 
regular estimator and at a wider range of 0. This suggests that the robust estimator has a 


smaller dependence on the emitter’s bearing compared to the regular weighted least- 
squares estimate. This also validates the effectiveness of the rejection processes in 
increasing the initial estimate’s accuracy in the presence of measurement outliers. From 
Figure 23, we see that the refined estimator’s performance is dependent on the initial 
estimator’s performance; therefore, any improvement to the initial estimate’s accuracy 


benefits the overall localization performance. 
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Figure 33. Root mean-square error ¢,,, of the weighted least-squares estimator 
(circle) and robust weighted least-squares estimator (diamond) versus @,. 
Results based on 3,000 trials with o;, =0 dB, N, =5, r,=200 m, and 

r,=2000 m. In both cases there were 2 random outliers present ata RDOA 


noise variance five times o;,,. 
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D. PERFORMANCE OF THE ROBUST LOCALIZATION IN THE 
PRESENCES OF MEASUREMENT OUTLIERS AND UNIFORM 
POSITION ERRORS 


In this section, we examine at the performance of the proposed robust weighted 
least-squares localization in the presence of measurement outliers and uniform position 
errors. To assess the performance of the localization with position errors we evaluated 
two Monte Carlo simulations, with and without measurement outliers present. Each 


simulation uses a 15 node ASLA configured network in which the emitter is located at 
o, =[2000 m, 30 deg]', r,=200 m, Tin =0 dB,, and all data points are the result of 
10,000 trials. 


A comparison of the robust and regular weighted least-squares estimator’s ¢,,. as 


a function of maximum uniform position error 0, is shown in Figure 34. In this case, 


there are no measurement outliers present. From the results, we see that the robust 


estimator does not improve performance over the regular estimator until a0, value of 0.4 


meters. Overall, this would suggest that the outlier rejection process is not beneficial in 


the small position error case given that there are no outliers present. 
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Figure 34. Root mean-square error ¢,,. of the weighted least-squares estimate (red) 
and robust weighted least-squares estimate (blue) versus é, . For 
comparison the no error case is shown in black. Results based on 10,000 
trials with oz, =0 dB, N,=5, r,=200 m, and w, =[2000 m, 30 deg]’. 
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Now, with measurement outliers present, a comparison of the robust and regular 
weighted least-squares estimator’s ¢,,. is shown in Figure 35. In contrast to Figure 34, 
the robust estimator shows a semi-fixed amount improvement over the regular estimator 
at all values of position errors simulated. This result suggests that the robust estimator is 
effective against measurement outliers but only marginally effective in the presence of 


positon errors. 
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Figure 35. Root mean-square error ¢,,. of the weighted least-squares estimate (red) 
and robust weighted least-squares estimate (blue) versus 6,. For 


comparison the no error case is shown in black. Results based on 10,000 
trials with o;,, =0 dB, N,=5, r,=200 m, and w, =[2000 m, 30 deg]. 


In this chapter, we developed a localization technique that is capable of 
approaching the CRLB. This technique consists of two sequential location estimators, a 
weighted least-squares estimator followed by a maximum-likelihood estimator to refine 
the initial estimate. We demonstrated that this technique is moderately susceptible to 
sensor position errors. Next, we developed a measurement outlier rejection process to 
compensate for these position errors as well as measurement outliers. The simulation 


results showed that the outlier rejection technique is effective against measurement 
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outliers and large position errors. In the next chapter, we provide a detailed description 


of the proposed collaborative signal collection technique. 
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V. ROBUST SIGNAL COLLECTION FROM AN AIRBORNE 
SEMI-STATIONARY NETWORK 


In this chapter, we develop the proposed collaborative signal collection scheme. 
This begins with an analysis of uniform and Gaussian sensor position errors and their 
effects on the ASLA formation’s array factor and main beam gain. Through this analysis, 
we are able to derive a novel signal estimator that uses sampled data and knowledge of 
the position error statistics to compensate for the sensor position errors. We further 
improve the collection’s robustness to interfering signals by adapting virtual filling and 
array tapering techniques to the proposed scheme. Simulations are used to support the 
developed theory and analyze its performance. 


A. EFFECTS OF UNIFORM POSITION ERROR ON THE ARRAY FACTOR 


In Chapter IV, we showed that uniform position errors negatively affect the 
localization performance. In this section, we discuss the effects of these errors on 


collaborative beamforming. 


From [63], we rewrite (28) to represent the array factor of a planar array as 
App (0.,,0,) > ae wel ? (85) 


where @, is the signal bearing,@, is the array steering angle, i.e., the direction of the 


main beam gain, @, is the i" node’s signal phase expressed as 


a, = B(x,sinO, + y,cos@,), (86) 
w, is the i" node’s complex weight expressed as 


W, a exp(—jA(xsin8,, + y,cos6,,)), (87) 
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M is the number of sensor nodes, and £ =2z/42. Incorporating uniform position errors 


o, and 6, (see (53)) into (86) yields 


G, = B((x,+5,,)sind, +(y, +6,,)cos@ ). (88) 


Separating the sensor position from its position errors, we get 


A, (8,,,9,,5,,5,) = ye wexp (ja, )exp(ja,, ), (89) 

where 
a,, = B(5,,sin8, + 5, cos, ) (90) 
is the i” node’s phase perturbation due to position errors. To steer the array to the 


intended signal’s direction @, we set @,, equal to @. Asa result, weexp( ja,;) =1, and 


the array factor becomes the expression for the main beam gain 


< y ei Ples,sind +5, s0084) (91) 


i=l 


Gye (9.6.0, ) = Ag (uF, ) 





6,4 =9, = constant 


where we assume that @ is known and treat it as a constant. By letting 


JB, Sin8, +6, ;cos8, 
Ne; =e ( e086, , we can express the mean of the phase error per node as 


1, = Elna} =], Pe eld, (5) f, (5) 4545, (92) 


By substituting f, (6, =F; iy (6,) = (26, ) into (92), we obtain 
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in(S, Bsind, )sin(d 6 
mpi oj EA) wy 
Dp t t 





and 


angle { My, } =0. (94) 
The variance of 77, ; 1s then obtained from 
o,,, =E\n.:\— i, (95) 
where the second moment E {72;} is expressed as 


. I f( 65, ;sin9, +6, ;c0s6, ) 


“fa, (8) fs, (5,) 46,45, 


Elm}= fe 








2 
1 5, po, jp(6,sin6,+6, ,cos0,) —j (6d, sind, +6, ,cos8, 
= | pP | pP e ( xi Tey a (Sei tT Cy 6. i 
26 ~5, 6, ts m 


(96) 
1 2 
op op 
=A 
Using (96), we rewrite (95) as 
o,. = I ae ae . (97) 


Assuming that 0, ; and 6, , are independent and identically distributed random variables, 


we have that 7,,; for i=1,...,M are also independent and identically distributed. Given 
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that Gy, (6,,6,,6,) is the sum of these M random variables, we express its mean and 


variance as 


M 
Gus = >on, =M pL, (98) 
i=l 
and 
M 
Conn = Do, =M (1-47, ) : (99) 


1. Effects of Uniform Position Error in the Special Case when the ASLA 
Formation Is Normal to the Target Source 


In the special case where the ASLA array has been reoriented (see Figure 18 (b)), 
the target emitter is normal to the array making 6 =@,,=Odeg. In this scenario, the 


expressions for the mean of the main beam gain is simplified such that (91) is now only a 


function of 6, and is written as 


Galo jay ee: (100) 


By letting 77, , = e’” we can express the mean of the phase error per node as 


5, 


ft, =E{fies} =| fe fs, (5,) 45,5. (101) 


By substituting f, (6,)=1/(26, ) into (101), we obtain 
Ay ry P 


magnitude { ji, } = sntef) (102) 
P 
angle {/7,|=0. (103) 


Along the lines of (98) and (99), we express mean and variance of the main beam gain as 
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Hg. =)4,, =M fi, (104) 
and 
M 
5, = 6, =M(1- A; ). (105) 
i=l 
2. Simulation Results on the Effects of Uniform Position Errors on the 


Array Factor 
To support these findings, three simulations were conducted using an N, =5 
ASLA configured network with a total of M=15 nodes each with uniform position 


errors. For these simulations, 2=1, the source emitter is located at a bearing of 


6 =45deg, r, = 200 m, and all data points are the result of 1,000 trials. 


A comparison between the simulated and theoretical values of £4, and Oy, as a 
function of 6, are shown in Figure 36 and Figure 37. From the results, we see that the 
simulated values closely adhere to the theoretical values, thus validating the derived 


expressions. 
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Figure 36. Mean normalized main beam gain of an N, =5 ASLA versus 6,. 
Results based on 1,000 trials withr, = 200 m, and 6. =45 deg. 
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Figure 37. Variance of the normalized main beam gain for N, =5 ASLA versus 
6,. Results based on 1,000 trials with, = 200 m, and 6 =45 deg. 


We now evaluate the effects of uniform position errors on the array factor. The 


average normalized array factor error y, (refer to (61)) as a function of maximum 
uniform position error 0, is shown in Figure 38. In this simulation, V,=5, each data 
point is the result of 1,000 trials, 2 =1, 7, =200m, and 6 =45 deg. From the results, 


we can see that the average y, rapidly increases with 0, , which means that even minor 


position(6, <l m) errors can have a significant effect on the array factor. This suggests 


that position errors of this nature if not compensated for will negate any benefits that 
collaborative beamforming can provide. This result is a key motivator for this research, 
as it is the objective of this dissertation to deliver an effective signal collection scheme 


(see Chapter V section C) in the presence of such errors. 
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Figure 38. Average y, for an N, =5 ASLA versus 6,. Results based on 1,000 
trials with 2=1, 7,=200m, and 8 =45 deg. 


B. EFFECTS OF GAUSSIAN POSITION ERROR ON THE ARRAY 
FACTOR 


In the previous section, we evaluated the effects of uniform position error on the 
array factor and the main beam gain of an ASLA configuration. In this section, we 


perform a similar analysis with Gaussian position errors. 


Similar to the uniform position error case, we wish to obtain an expression for the 
mean and variance of the main beam gain; however, using the same methodology 
becomes intractable for the Gaussian random variables. As an alternative, we instead 
evaluate the main beam gain after the ASLA formation has been reoriented as shown in 
Figure 18 (b). In this configuration, we assume that 9 =@,,=0 deg. This assumption 
allows us to study the effects of Gaussian position errors independent of errors in the 


bearing estimate 6 . As aresult of this assumption, the main beam gain from (91) can be 


expressed as 


Tid 


Glo => es (106) 
where €”” is known as the wrapped Gaussian random variable [84]. By letting 


gi , we can interpret the main beam gain as the sum of M wrapped Gaussian 


7 . " 2 
random variables , each with a mean of zero and variance o; = ( Bo,) [84]. The 


probability density function of the wrapped Gaussian random variable is given by [84] 





2 
: 2nk 
fur, (Wes) = — Dope °XP ve) (107) 
& g& 


The mean values of the magnitude and phase, respectively, are given by [84] 


magnitude { ,,} = gical (108) 
angle{,,=0. (109) 

The variance is given as 
sere), (110) 


Since the main beam gain is the sum of M terms, similar to (98) and (99) we express the 


mean value of its magnitude and phase as 


magnitude { 444, | = Me (7) (111) 


angle { 14,,, | = 0 (112) 


with a corresponding variance of 


o?,=M [1 agi)? } (113) 
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1. Simulation Results of the Effects of Gaussian Position Errors on the 
Array Factor 


To support the theory presented in the previous section, three simulations were 


conducted using an N, =30 ASLA configured network, with Gaussian position errors. 
In these simulations, the source emitter is located at a bearing of 6 =Odeg, 2 =1, 


r, = 200 m, and all data points are the result of 1,000 trials. 
A plot of the simulation and theoretical normalized values of magnitude { dig) as 


a function of oO; =/( Bo, y is shown in Figure 39. From the plot, we see that the 


simulation results closely follow the theoretical. The results indicate that the main beam 


gain monotonically decreases as oO. increases. Similar to the uniform position error case, 


we see that the errors significantly reduce the effectiveness of collaborative 


beamforming. 
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Figure 39. Normalized values of magnitude { /1,,,} versus position error variance 


ee (Bo, .: . In this simulation V, =30, 6 =Odeg, 2=1, r,=200m, 
and all data points are the result of 1,000 trials. 
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Similar to the previous simulation, a comparison between the simulation and 

: 2: 5. se 
theoretical values of angle{,,,} as a function of o; = ( Bo,) is shown in Figure 40. 
From the results, we see that the simulation results closely follow the theoretical value, 
with the angle { ne values approximately equal to zero at all values of oe simulated. 


This suggests that the effects of Gaussian position errors have no effect on the mean 


value of the main beam phase angle. Since position errors do not affect the angle{ ais ‘ 


it will have a prominent role in the signal estimation described in the following section. 
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Figure 40. Values ofangle { /,,,} versus position error variance * P! | In this 
simulation V, =30, 8 =Odeg, 2=1, r, =200 m, and all data points are 
the result of 1,000 trials. 


c. SIGNAL COLLECTION FROM AIRBORNE SYMMETRIC LINE ARRAY 
NETWORK IN THE PRESENCE OF UNIFORM POSITION ERRORS 


In the previous section, we showed that each node’s phase perturbations due to 
position errors have a mean value of zero. Given that the ASLA configuration contains 
three stacks of N, sensor nodes each, we consider a snapshot of the complex signal 


Ry 


received by each node in a particular stack, which we model as the true signal received at 
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that stack with phase perturbations due to position errors. By using the mean value of the 
in-phase and quadrature components of all the measurements from a particular stack, we 
can estimate the true signal magnitude and phase at that stack. Since position errors do 
not affect the mean value of the phase estimate, this method will recover the signal 
magnitude and phase despite the presence of position errors. Given this, we derive an 
expression for the mean of the in-phase and quadrature components at each stack. This 
expression is used to estimate the phase and magnitude based on the sample mean of the 


received complex snapshots. 


Considering an ASLA configured network with uniform position errors, we 


represent a snapshot of the complex signal received at the i" node of the center stack as 


id, JB 6, ;8inO, +6, ;cosd, 
1o, = Viele ), (114) 


where i=1,...,N,, and V, and @ are the signal’s magnitude and phase, respectively. 


Using the same assumptions used to derive (106), we express (114) as 


Voy = Veoh) (115) 


h 


where p,; = 80,, is a uniform random variable representing the i’ node’s phase 


perturbation due to position errors in the y-direction. The resulting complex signal is 


then represented as two random variables [85] 


Vee gg eI ss (116) 
where vy. . is the i” node’s in-phase component expressed as 
I 6, il 
¥.,4 =V,cos(g, + 2.) (117) 


and Piss is the i” node’s quadrature component expressed as 
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¥.,i =Vsin(g, + p,,)- (118) 


We consider y, and 7, a8 two random variables representing a given stack’s received 
in-phase and quadrature components, respectively. The random variables y, ; and Vict 
can be interpreted as the i" instance of the true y., and Ve, that are perturbed by node 


position errors. Also, the mean values of all the instances would then be the sample 


means ofy, and Fox Since we have shown that position errors do not affect the mean 
value of the signal phase, we can use the sample mean of y, and ¥,, to derive an 
estimate of the true signal phase. Overall, by deriving an expression for the mean of y, 


and Yo.» we can also determine a way to estimate the true signal phase. 


With a goal of derive an expression for the mean of ve and y, for a given stack, 
we begin by deriving the probability density function of Y,.. By substituting p, =% +p, 


in (116), we get 
Yeo =Vsin(p,). (119) 
The probability density function of 9,can be expressed as 


1 
A (0) =F 55” $, - BS, < p, <4 + P6,. (120) 


We next substitute «, = sin(p,) in (119) to get 
¥,. =V,K,. (121) 


The probability density function of the random variable «, is expressed as 
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1 1/(2B6 
fe Te /dp fo. (P,) eA (122) 
s s p.=8(k,) 1-(x,) 
where 
g'(x,)=sin' Ky» (123) 
and 
ldx, /dp, Sees cos( sin” K,) =J1-(x,), (124) 


for sin(¢, - B6,)< K,< sin(¢, + B6,). Finally, the probability density function of the 
random variable 7, =V,x, is expressed as 

Ky 

V, 


Sr, (7, )=> he 


t 





(125) 
1 1 
= , Vsin(¢g-Bo,)<y7-. <V,sin(¢ + Bo, ). 
5B, Wore, (9, - 85, )<rc, <Vsin(9, + 86, ) 
t o 
We then account for measurement noise using the expression 
d. LG, +L yen (126) 
and 
C= ones Quen > (127) 


where J yyy and Qjycy are independent Gaussian random variables, representing the 
complex additive white Gaussian noise. Using these expressions, we derive an estimate 


of V, and ¢ from the mean of all the center stack’s complex in-phase and quadrature 
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measurements. The mean of each stack’s in-phase and quadrature components can be 


expressed as 
Ho, =E{I}=El7¢ }+E {Loven} (128) 


He, = E{Q,$= E{¥0,}+ E {Oana} - (129) 


With regard to the quadrature component, sinceQ,,,. is zero mean, Hc, can be re- 


written as 


He, =E}7c,| 


ie 1 ior 


a 
«2a, epee ISSO 


Y dl x 
[es ~ 2 dy,+] 5 = 2 Chey (130) 
Cav, ae on sd NLA V5 ; 


a (a 
5 <|als>te 








p 


ibe 1 05 


a 
via, rn, SP a% 


where @, = fo,, Ay =V,sin (4, -a,) , and By = Vsin (4, +a,) . We define the mean of 


the in-phase component similarly as 
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qr dY¢,; peal 
a pa ie eee eae, ve 





where A, = V.cos(¢, -a,) and B, =V,cos (4, + a,). By solving (130) for ¢, we obtain 


an expression for the true signal phase that is a function of He,» V, and @,. Then using 


A 


the sample mean He, in place of the true mean Hc,» We obtain the phase estimate ¢,. 


The derivation of the phase estimate ¢, for the case of \4| <a/2-a, is shown as 


A. exe Bo 1 Yc 
~ I, 2a, I? 7 i dV c,1 


a Bo 


=| WV? -(r i (132) 


PL Ao 


=a Wea) Wee (| 











By squaring both sides and expanding Ag and Bg terms, we get 
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(20,f4c,) -[ -(Vsin( (¢, a, ie (Vsin( (4, £Q, yy. (133) 


Expanding the right hand side and applying the trigonometric identity 


cos’ @-sin* 0=cos20, we get 


(20,Ac, ) =) (1 +cos(2¢,)cos (2a, )) 


(134) 
—2,/V,*cos* (¢, — a1, )cos” (4, + a,) 
By taking V;* out of the square root and bringing V,’ to the left hand side, we get 
2 
2a fi.) =V,(1+cos(2¢,)cos(2a 
( P i) ( ( ) ( »)) (135) 


—2V? 





cos( 4, —a,,)cos(¢, +a,) 





By applying trigonometric identity cos xcos y =(1/2)[cos(x+y)+cos(x—y)] to the 


absolute value terms, we get 











2 A 2 
(7A =1+ cos(2¢,)cos(2a, )- \cos(2@,) + cos (2¢, )). (136) 
We perform the same derivation for the mean of the in-phase component to obtain 
2a ul ‘ 
(AeA. =| —cos(24,)cos(2«, )-|cos(2@, ) -cos(2¢, ) ; (137) 
Taking the ratio of (136) and (137) yields 
He, ? 1+-cos(24,)cos(2a, ) -|cos(2a, ) +cos(2¢, 
Sy SS SS 138 
He, 1—cos(2¢,) cos(2a, )- [eos (2c, )—cos (24) ee 
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2 
To solve for ¢, we substitute, 1, = ie | Be ,a= cos(2a, ) and b =cos(2¢,) in (138) 
to get 


Z 1+ab—|a+b| 


= —___—__. 139 
Bs 1-ab—|a—b| o27) 


By taking only the positive solutions for the absolute value terms and solving for b, we 


get 


pose niles, (140) 
1+ yu, 


Since b = cos(2¢,), we can solve (140) for the phase estimate 


g, =(1/2)cos™ (EH) (141) 


For the other two cases of (130), repeating the above steps yields the same expression for 
the phase estimate d. By substituting d, into (137) and solving for V;, we get the signal 


magnitude estimate 


























p= 
20 phe, b|<-a 
, | = 
\/1-cos(2¢, )cos (2c, )- cos(24,)—cos(2a, } Zin s* 
2et, fl 2» clgleZag 0 
cos(¢, -@, ||+[cos(¢, +a, Poe ‘ 
da i 
aK, lal> 5-2 


1—cos(2¢, )cos(2a, )- cos(2¢, )—cos(2a, ) 
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With each stack in the array using the same estimation technique, we express their 


complex signal estimates as 


A A 


© =V,,exp( jb), 2=12,3, (143) 


where V,, and d,, are the left, right and center stack’s magnitude and phase estimates, 


respectively. Each stack’s Q, is then treated as the complex output of a single node in a 


phase shift beamforming array. The final main beam response is then expressed as 


3 
Oey Os (144) 
n=l 
1. Signal Collection in the Presence of Thermal Noise and Uniform 


Sensor Position Error 


To demonstrate the performance of our signal collection and estimation scheme, 
we present the results of five Monte Carlo simulations in Figure 41 through Figure 45. 


The results of these simulations focus on the root mean-square error for signal magnitude 
(é,) and phase ce) estimate. For all five simulations,2 =1m, the source emitter is 
located at w, =[2000 m, 45 deg]', r,=200 meters, and all data points are the result of 


10,000 trials. Note a key parameter in the following simulations is signal to noise ratio 


(SNR) in units of dB. This quantity is a measure of the ratio between signal power P. 


and noise power P, and is expressed as SNR = 10log,,(P./P,). 

A comparison of the scheme’s ¢, and ¢, as a function of per node SNR at 
different RDOA noise variance o;, levels are shown in Figure 41 and Figure 42, 
respectively. The target bearing estimate Gis derived using (75). From the results, we 


see that the effect of the per-node SNR on ¢, and ¢, is not as strong as that of the 6 


estimate. This suggests that the estimator is heavily dependent on the accuracy of the 


localization. 
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Figure 41. Root mean-square error ¢, at different levels of RDOA noise versus 
SNR per node. Results based on 10,000 trials with V, =10,7, = 200 m, 
5, =.05 m, and @, =[2000 m, 45 deg]’. 
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Figure 42. Root mean-square error ¢, at different levels of RDOA noise versus 
SNR per node. Results based on 10,000 trials with V, =10,7, = 200 m, 
5, =.05 m,and @, =[2000 m, 45 deg]’. 
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Two plots comparing the scheme’s ¢, as a function of maximum uniform 


position error 6, and the number of nodes per stack N, as a function of o;, are shown in 
Figure 43 and Figure 44, respectively. In Figure 43, we see that uniform position error 


doesn’t affect ¢, until values of 5,>0.1m in the low op, case(o@ =—10 dB). 


Conversely, it causes a near linear increase in Ce for all values of é, simulated. In 
Figure 44, the results show that the benefits of increasing N; diminishes when 
N, >10 nodes. Unlike the results from Figure 21, where increasing the number of nodes 


improved performance, here increasing the number of nodes per stack is less effective 


above a certain value. 
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Figure 43. Root mean-square error ¢, at different levels of RDOA noise under 
increase maxed uniform position error 0,,. Results based on 10,000 trials 


with, NV, =10, r, =200 m, SNR =5 dB, and @, =[2000 m, 45 deg]’. 
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Figure 44, Root mean-square error ¢, at different levels of RDOA noise under 
increase maxed uniform position error V,. Results based on 10,000 trials 
with, 5, =1/10 m,r, =200 m, SNR =5 dB, and @, =[2000 m, 45 deg]". 
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A plot of the proposed scheme’s ¢, as a function of per node SNR with and 
without signal estimation is shown in Figure 45. Each case was also evaluated at two 
different values of o;,,. The “with estimation” case represents our proposed scheme, i.e., 


collaborative beamforming in conjunction with the proposed signal estimator. The “no- 
estimation” case represents collaborative beamforming without the use of signal 


estimation. In this simulation, the source emitter is located at@, =[2000 m, 45 deg]', 
N,=10, A=1m, é, =1/10m, and all data points are the result of 10,000 trials. In the 


Onn =—15 dBease, we see that the benefits of the signal estimation are evident at all 
values of SNR simulated. As evident in the oj, =-5 dB case where the scheme 


underperforms the no-estimation case, we again see that the scheme is highly dependent 
on the accurarcy of the localization phase. This dependency is derived from the scheme’s 
assumption that the target signal’s bearing is perfectly normal to the array, i.e., the 


@. =0 deg, when in actuallity @ it is a small random quantity based on the accuracy of 


the localization. 
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Figure 45. Signal ¢, of proposed scheme compared with the no-estimation case 


versus SNR per node values. Results based on 10,000 trials with 
5, =1/10m,r, = 200 m, and @, =[2000 m, 45 deg]’. 
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An illustration of the mean value of the normalized gain of the array as a function 


of signal frequency for different values of 0, is shown in Figure 46. For comparison, we 
also include the no-estimation case. In this simulation, the source emitter is located at 
o, =[2000 m, 45 deg]", N,=10, A=1m, 56,=1/10m, og, =0 dB, and all data points 
are the result of 10,000 trials. As frequency increases, so does f£. From (115), we see 
that phase error due to 6, is proporitional to 6. As #=27f,/c, an increase in 


frequency will magnify the negative effects of the position errors. From the results, we 
see that allowing these position errors into beamforming (no-estimation case) results in a 
response similar to that in Figure 36 where the main beam gain is significantly reduced as 
position errors increased. When we incorporate the proposed signal estimator, an 
improvement over the no-estimation case is apparent for all simulated frequency values. 
This suggests that the inclusion of signal estimation in collaborative beamforming can 
effectively mitigate array phase pertabations due to the uniformly distributed positioin 


errors. 
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Figure 46. Average signal normalized gain at different levels of RDOA noise 
versus signal frequencies. Results based on 10,000 trials with a. = 0 dB, 


5, =1/10m,r, =200 m, and @ =[2000 m, 45 deg]". 
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D. SIGNAL COLLECTION FROM AIRBORNE SYMMETRIC LINE ARRAY 
NETWORK IN THE PRESENCE OF GAUSSIAN POSITION ERROR 


In the previous section, we developed a signal estimator to combat the effects of 
uniformly distributed position errors. In this section, we derive a similar estimator for the 


case of Gaussian distributed position errors. 


In a fashion similar to (116), a snapshot of the complex signal received at the i” 


node of the center stack is represented as 


No = Ne, tc (145) 


where 


Nc, =V,c0s(, +P, ) (146) 


with 9, = Bo, with é, being zero mean Gaussian position error with a variance of 07; 


(see Chapter III) and 


Nc, = JVsin(g, + P,). (147) 


Similar to the uniform position error case, we will derive an expression for the mean 
value of the in-phase and quadrature components for each stack and then use these 
expressions to derive an estimate of the signal phase. Also similar to the uniform random 


; : : sth - 
variable case, each node’s n,.. andn,. . ina stack can be interpreted as the i" instance of 
p) C,,i Co wl 


the in-phase random variable n,, and the quadrature random variable Nc, - 


Focusing on the in-phase component first, the expression @, + 2, can be expressed 


as a Gaussian random variable 0, with a mean of 


H, =§, (148) 


and variance of 
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2 
o,= (Pax). (149) 
By substituting 0, into n, , we get 
Nc, = cos(d, ). (150) 


To derive the resulting probability density function of 1, , we first define its cumulative 


distribution function as 


Fy, (ne, )=Pr[ 9(5,)<ne |= J tu, (6) 46. (151) 


a(d, )sne, 


where g (6,) = COs (6,) . Since a Gaussian random variable is define between —co and , 


6, can be defined by the inverse transform of g (6, ) over the interval [84] 


2nk + cos (Mg, ) £0, < 2m(1+k)—cos"' (ne, ) (152) 
The probability density function of 1, is then expressed as [81], [86] 


do 


vee 


dn 


Fig (Me, )= Line 








(fs (27k +cos | (n., )) + fu (27(1 +k)-cos | (n,, ))) , (153) 


C, 


where 


= (154) 
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By expanding (153), we get the final expression for the probability density function of 


Nc, as 








1 1 4-4 8-4 
fa (e)= Spe po mal Io }} (155) 


where 4, =2nk+cos' (n., ) and B, =2n(1+k)-cos' (n., ). 
To support the validity of this expression, a comparison of the histogram of 


Nc, =V,cos (6, ) with the response of (155) over the same interval is shown in Figure 47. 
For this simulation, 10,000 instances of 0, generated to form a histogram with 50 bins, 
o, =(22/ 5)’, and fy. (n., ) is calculated with k =0,1. As observed from the results, the 


simulated histogram closely follows the theoretical values of (155). 
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Figure 47. Comparison of measured and theoretical probability density function 
(pdf) of n, . Signal phase ¢, = 0 and measured pdf based on histogram 


using 10,000 points and 50 bins. 
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1. Mean Value of the In-Phase and Quadrature Components with 
Gaussian Position Errors 


With the probability density function of 1, defined as a function of g and om 
we can derive an estimate of ¢, at a given stack based on the sample mean of 7, and 
knowledge of om Since it is difficult to derive an expression for the mean using the 


probability density function of 1, , we take an alternative approach using Euler’s formula 
iby at 
e * =cosd, + jsino, . (156) 


We will first consider the case where 6 , is Gaussian with a mean value of zero and a 


variance of O.. Note that the mean value of 6, is indeed equal to @, (see (148)), but we 


will first examine this simpler zero mean case in order to derive the non-zero mean case. 


The mean of the cosine term is then given as [87] 


E{cosd,\ = exp(-o; /2), (157) 
and 


E{sind,} =0. (158) 
Using the trigonometric identity 
1 
cos” é, =5(1+c0s(26,)), (159) 
the variance of cosd, can be expressed as 
var {cosé, | = E{cos’d,—E{cosd,}". (160) 
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By substituting (157) and (159) into (160), we get 
1 -(20, ye -o; 1 -o; 2 
var{oosd,}=>(I+¢ -e =5(I-e ] : (161) 
Similarly, the variance of sind, is given by 
. + 2 . 2 1 207 
var {sind, | = E {sin 6,|\-E{sino, \ =5(I-e \. (162) 


Finally to consider the case where 6, is not zero mean, where 6, is a Gaussian random 
. . = . 2 . . = _ 
variable with mean wv, = ¢, and variance equal to o,. If this is true, then 0, =0, —H, 


where 6, is the zero mean case of 6,. Thus, using the trigonometric identity 


cos(x + y) =cosxcos y—sinxsin y , the mean of cos(6, } can be expressed as 


E{cos 6, = Eteos( 4, + Ms, ) = E{cos OpCOSL, —sind,siny, | = cosy, (163) 
Similarly, 
E{sin 5,}=E{sin(5, + u,)}=e 7 ?sin es (164) 
Sincen, = cos(d, ) , the mean of 7, is given as 
E{n,,} =e" cospt,. (165) 
Following the same derivation, the mean of Nc, is expressed as 


E\ng, } =e "sing, (166) 


97 


2. Phase Estimates Using Sampled Base Band Signals 


Given that the signal phase is perturbed by sensor position errors, as shown in 
Chapter V section B, the mean value of the signal’s phase remains largely unaffected. 
This makes the sample mean of a given stack’s signal phase a suitable parameter from 


which to estimate the true signal phase. 


One means of estimating the signal phase is by solving (163) or (165) for @, to 


obtain the phase estimate for the center stack 
scos'(& /e*?)=sin'[f se), (167) 
tc H,., Hig, 


where i, and ie are the signal’s in-phase and quadrature component’s sample mean 
1 Q 


for the center stack. This method, although valid, only uses the samples from either 
component. The preferred approach is to use both base band signal samples by taking the 


ratio of (166) and (165) as 


Ejn, e7 "sinus, singly 
a} ns (168) 
E{n,} ¢% Cosp, —COSHs, 





Since UW; =@., an estimate of ¢ . is expressed as 
6, t,C t,c Pp 
i Bis 
,, = tan [Se (169) 


This same phase estimation method is then used for the other two stacks in the array. The 
final signal phase estimate d, is then the average of these three estimates. Since the 


ASLA array is symmetrical about the y-axis, taking the average of the estimates will also 
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compensate for the minor errors in the arrays reorientation, i.e., errors in the estimate @. 


The estimated complex signal for the n™ stack is then expressed as 


©, =exp( ig), (170) 
We then sum the estimates of all three stacks to obtain the final main beam response as 


fe) “50; (171) 


3: Performance of the Signal Collection Scheme in the Presence of 
Gaussian Sensor Position Errors 


To evaluate the proposed scheme’s performance in the presence of Gaussian 


sensor position errors, we conducted a series of Monte Carlo simulations. The results of 
these simulations focus on the root mean-square error for signal magnitude (é,) and 
phase (é,) estimate. For comparison, the no-estimation case, i.e., beamforming without 


signal estimation, is also shown. For each of these simulations,/=1m, the source 
emitter is located atw, =[2000 m, 15 deg]", r,=200 meters, andN,=10. All data 


points are the result of 10,000 trials. 


A comparison of the scheme’s root mean-square error of normalized magnitude 


¢, with and without signal estimation as a function of the standard deviation of Gaussian 
position erroro, is shown in Figure 48. In this simulation, to isolate the effect position 
error, the target signal’s @ is known a priori. From the results, we can see that the 
estimated case outperforms the no-estimation case at all values of o, simulated. This 


suggests that signal estimation provides an effective means of reducing the effects of 


Gaussian position errors. 
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Figure 48. Root mean-square error of normalize magnitude versus the standard 
deviation of Gaussian position error o,,. Results based on 10,000 trials with 


N,=10, r, =200 m, and , =[2000 m, 15 deg]". 


The scheme’s ¢, as a function of Gaussian positon error o, is shown in Figure 
49. Similar to the previous simulation, we assumed the target signal’s @ is known a 
priori. From the results, we observe no noticeable difference in the Ce between the 


proposed signal estimation case and the no-estimation case. Similar to the uniform 
position error case, this suggests that the proposed signal estimation technique does not 
improve the accuracy of the main beam signal phase, but as seen in Figure 48 does 


improve its magnitude gain. 
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Figure 49. Root mean-square error phase ¢, versus the standard deviation of 
Gaussian position error 0,. Results based on 10,000 trials with V, =10, 


r, =200 m, and @, =[2000 m, 15 deg]’. 


The scheme’s &, as a function of RDOA noise variance 0%, =c’o; is shown in 
Figure 50. In this simulation, the scheme is evaluated at two levels of position error with 


0, =0and o, =1/5m. From the results, we can see that the proposed scheme and no- 


estimation case have the same performance when no position errors are present. When 
position errors are introduced, we see that the proposed scheme outperform the no- 
estimation case. Of interest, we see that the scheme also outperforms the same case with 


no position error for values of RDOA noise greater than 2 dB. This would suggest that at 


‘ 9 ge , a 
a given 07, level, position errors become less of a factor than errors in @.. 
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Figure 50. Root mean-square error of normalize magnitude versus RDOA noise. 
Results based on 10,000 trials with, V, =10, 7, =200 m, and 


@, =[2000 m, 15 deg]". 


A comparison of the scheme’s Ce as a function of Ope is shown in Figure 51. 


Similar to the previous simulation, the performance was evaluated at two levels of 


position error, 0, =0 and o, =1/5 m. Similar to Figure 49, we can see that the proposed 


scheme and no-estimation case have the same performance with or without position 
errors. Also, similar to the previous simulation, we see that the cases with position error 


beat the no position error cases at values of RDOA noise greater than 2 dB. 
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Figure 51. Root mean-square error Phase ¢, versus RDOA noise. Results based on 


10,000 trials with, N, =10, r, =200 m, and @ =[2000 m, 15 deg]. 


E. GRATING LOBE SUPPRESSION VIA VIRTUAL FILLING 


From (55) and (56), we know that increasing the distance r, between the outer 


and center stacks (see Figure 6) will provide increased localization performance, but 
increasing this distance essentially increases the inter-node distance d in (28). If the 
inter-node distance of an array is greater than 1/2, multiple grating lobes will appear in 
the array factor [71]. Since grating lobes have the same magnitude response as the main 


beam, they will significantly reduce the array’s interference rejection capabilities. 


To combat the effects of grating lobes, we propose the use of virtual filling. 
Given that the directions of arrival of all signals are known, virtual filling interpolates the 
response of a virtual uniform linear array based on signal estimates from the real array 
[76]. In virtual filling, any gap between array nodes larger than 1/2 is filled in with 
virtual complex signal data. 
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With the ASLA symmetrically populated only along the x-axis and centered on 


the origin, the position of the n™ virtual node is expressed as 


1 gare Ie 1, = 0,1 «2,N M1 (172) 
where the virtual inter-node distance d, =2/2 and N, is the total number of virtual 


elements given as 


N,=—*41. (173) 


To fill in the complex data of the virtual array, we begin by populating the vector 
A, from (51) with the signal estimates Q,. Using the least-squares method described in 


Chapter II, an estimate of each signal’s magnitude and phase is obtained. The response 


of the 7" virtual node is then expressed as 


0,,=387,e%eiMesint 74) 


vn t,i 


h 


where N, is the number of signals present,i=1,...,N,, and V,, ,¢,, @, are the i' 


signal’s magnitude, phase, and bearing, respectively. Once filled out, the ASLA then 
becomes a virtual uniform linear array. Since the internode distance is now< 1/2, the 
presence of grating lobes is significantly reduced. Furthermore, various tapering methods 
(see Chapter II) can be applied to manage the remaining sidelobes. 

With no errors present, the array factor of an N, =5 ASLA configured network is 
shown in Figure 52. For this simulation, 2 =1m and r,=2002. From the plot, we can 
see that with r,>>A/2, there are numerous grating lobes and side lobes in the array 


factor. 
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Normalized Power (dB) 





Figure 52. Normalized Array factor of an NV, =5 ASLA phase shift beamformer 
before virtual filling process. Simulation based on r, = 200 m, and 
o, =[2000 m, 45 deg]. 


Now, with virtual filling and tapering applied, the array factor of the ASLA is 
shown in Figure 53. With the virtually filled array having no inter-node spacing greater 
than 2/2, we can see from the plot that all the grating lobes have been eliminated. By 
applying a binomial taper (see Chapter II) the sidelobes have also been suppressed, 


resulting in a strong main beam lobe. 


-30 30 


60 / \ 60 
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Figure 53. Normalized Array factor of an NV, =5 ASLA phase shift beamformer 
with virtual filling process with binomial tapering. Simulation based on, 
r, =200 m, and o, =[2000 m, 45 deg]. 
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A comparison of the signal estimation scheme with and without virtual filling is 
shown in Figure 54 and Figure 55. In this simulation, there is an interfering signal in 


addition to the target signal, 4=1m, uniform position error 6, =1/20 m, and the ASLA 
formation consists of N,=10 with r,=100m. The target signal’s 6 =0 deg, V, =1, 


and ¢, =60 deg. The interfering signal’s 9 =55 deg, V,=1, and @ =25 deg. 


Ss 


The magnitude and phase estimate values of the main beam for 100 trials without 
virtual filling are shown in Figure 54. From the results, we can see that the interfering 
signal corrupts the target signal, thus resulting in an erroneous main beam response. In 
contrast, the main beam signal values using virtual filling are shown in Figure 55. From 
the results, it is clear that the virtual filling has reduced the effects of the interfering 


signal and has improved the signal estimate. 
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b) Phase estimate ‘without virtual filling 
Figure 54. Signal Estimation without virtual filling. Simulation based on, 
r, =100 m, target signal’s 0 =0 deg, V. =1, and ¢ =60 deg. The 
interfering signal’s 6 =55 deg, V,=1, and ¢ =25 deg. The true signal 


phase and magnitude is shown in red. 
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Figure 55. Signal estimation with virtual filling. Simulation based on, 7, =100 m, 
target signal’s 6 =0 deg, V.=1, and ¢ =60 deg . The interfering signal’s 
6, =55 deg, V,=1, and ¢, =25 deg. The true signal phase and magnitude 

is shown in red. 


A comparison of the main beam’s angular root mean squared error ¢, for the 
proposed scheme with and without virtual filling as a function of maximum uniform 
position error é, is shown in Figure 56. In this simulation, there is one signal-of-interest 
and one interfering signal, the target signal’s 8 =O deg, V,=1, and ¢ =60 deg, the 
interfering signal’s 6. =55 deg, V,=1, g =25 deg, and 7, =100 m; and all data points 


are the result of 1,000 trials. From the results, we see in the no virtual filling case that the 
benefits of the signal estimation are nullified by the presence of an interfering signal. In 
contrast, when virtual filling is applied, the interfering signal is significantly reduced, 
thus improving the signal phase estimate. Overall, this result demonstrates the proposed 


scheme’s effectiveness in the presence of both position errors and an interfering signal. 
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Figure 56. Proposed scheme’s ¢, with virtual filling (blue) and without (black) 
versus maximum uniform position error 6, in the presence of one 
interfering signal. In this simulation, 7, =100 m, the target signal’s 
0 =0 deg, V,=1, and ¢ =60 deg, the interfering signal’s 9, =55 deg, 
V,=1, and ¢@ =25 deg, and all data points are the result of 1,000 trials. 


In this chapter, we studied collaborative signal collection in the presence of noise, 
sensor position errors, and interfering signals. We developed a signal estimation method 
that utilizes sampled array data and knowledge of the positional error statistics to combat 
these errors. We further enhanced its robustness against interfering signals by using a 
virtual filling and tapering techniques to minimize the array’s grating lobes and side lobe 
response. Using simulation, we were able to validate the proposed scheme’s 
effectiveness in the presence of such errors and achieve the intended objective of this 


dissertation. 
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VI. CONCLUSION 


Performing collaborative beamforming from an elevated, mobile wireless sensor 
network requires a coordinated interplay of many different signal processes and 
technologies. The objective of this research was to maximize signal collection 
performance in the presence of various signal and sensor related errors from an elevated, 
mobile wireless sensor network. To accomplish this objective, we proposed a signal 
collection scheme that exploits an elevated, mobile network to maximize the 


collaborative collection of a target signal. 


The proposed scheme begins with an accurate source localization technique, 
which is used to determine the signal’s position and bearing. This information was then 
used by a collaborative beamformer to maximize the collection of that target signal. Both 
of these phases are enhanced by the versatile nature of an elevated, mobile network. The 
mobility allows for the reconfiguration of the sensor network’s topology to help create an 
ideal sensor-target geometry. This geometry allows for optimal localization estimates to 
be obtained. The ability of an elevated network capable of creating sensor stacks allows 
for robust signal collection operations in the presence of uniform and Gaussian sensor 


position errors. 


In the localization technique, to obtain precise localization in the presence of 
noise and uniform position errors, we proposed the use of two sequential location 
estimators. This technique consists of an initial weighted least-squares estimate followed 
by a maximum-likelihood estimate. The second estimate uses the initial estimate to 
reorient the network. This reorientation creates an ideal sensor to target geometry [17]. 
This new orientation minimizes the geometric dilution of precision [52], minimizing the 
maximum-likelihood estimator’s error variance. Simulation results showed the proposed 
localization technique to be efficient, i.e., the variance of the estimation error approaches 


the Cramer-Rao lower bound. 


To enhance the localization technique’s performance, we developed a 


measurement outlier rejection process that mitigates the effects of measurement and 
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sensor position errors. This technique uses a combination of single case diagnostics and 
the Mahalanobis distance to identify and remove erroneous measurements from the least- 
squares estimate. Through simulation, we demonstrated this technique to be effective in 
the presence of both measurement and position errors. For the case where only position 
errors are present, the technique was found to be effective only for high levels of position 


errors. 


For the signal collection technique, we began by analyzing the effects of sensor 
position errors on the array factor and its main beam gain. We then derived an 
expression for the effects of both uniform and Gaussian position errors on the array’s 
main beam gain. Through simulation, we validated our results and showed that the mean 
value of the main beam signal phase was unaffected by position errors. To enhance the 
signal collection performance in the presents of these errors, we developed a signal 
estimator for both uniform and Gaussian position errors. Simulation results yielded an 
improvement in array gain of approximately 37 percent for the standard deviation of 


position error values greater than 0.4 m. 


Taking advantage of the concept of a unique elevated, mobile wireless sensor 
network realized through the use of multirotor UAVs, our scheme used two existing 
localization techniques to deliver a precise source location estimate. Using this 
information with statistical knowledge of the sensor position errors in a signal estimator, 
we derived a novel collaborative signal collection scheme. This scheme was shown to be 
capable of collecting and amplifying a target signal in the presence of such errors. With 
all these techniques in concert, the objective of this dissertation to maximize signal 


collection in the presence of various signal- and sensor-related errors was achieved. 


A. SIGNIFICANT CONTRIBUTIONS 


With the objective to maximize signal collection performance in the presence of 
various signal and sensor related errors, a novel scheme for robust signal collection from 


an elevated, mobile network was proposed. 


To enable the signal collection, a localization step was required. To accomplish 


this task, a two-stage localization technique was proposed. Using an optimal ASLA 
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sensor formation for target localization, we minimized the network’s geometric dilution 
of precision, effectively lowering the estimate’s Cramer-Rao lower bound. By utilizing 
two unique location estimators in sequence, we created a hyperbolic location estimate 


capable of approaching the Cramer-Rao lower bound. 


To enhance the localization robustness, a measurement outlier rejection process 
was proposed. Using the Mahalanobis distance as a measure of solution divergence, we 
were able to identify and remove erroneous measurements from the location estimate. As 
a result, the localization estimation technique demonstrated increased resilience to both 


measurement and sensor position errors. 


To better understand the effects of position errors in collaborative beamforming, 
an expression for uniform and Gaussian sensor position error effects on an ASLA 
formation’s main beam gain were developed. The main result was the validation that 
position errors in the ASLA formation have no effect on the mean value of the signal 
phase, which was instrumental in the development of a novel signal estimator used to 


compensate for such errors. 


Finally, a signal collection scheme that combines a novel signal estimator and 
collaborative beamforming to successfully collect a target signal despite the effects of 
sensor position errors was proposed. Compensating for both uniform and Gaussian 
errors, the signal estimator uses sampled array data, sensor stack formations, and 
knowledge of the positional error statistics to accomplish this task. Through simulation 
we were able to demonstrated significant improvement over the conventional 
beamforming. To further enhance the collection capability for noise and interference 
rejection, we also adapted the use of virtual filling and tapering to manage the array’s 


grating lobes and sidelobes. 


B. FUTURE RESEARCH 


With the research and development of localization and signal collection 
consistently changing and growing, there are a number of possible extensions to the 


research presented in this dissertation. 


111 


The key error parameter evaluated in this research was sensor position errors 
modeled as a uniform or Gaussian random variable. The assumption of specific 
probability density functions affects many aspects of each of the techniques and methods 
developed in this research. Further research focused on a general approach to treat 
position errors independent of a probability distribution should be considered in order to 


fully investigate their effects on signal collection. 


The robust source localization technique presented here focused on a single case 
diagnostic. This method assumes that single data points act alone to offset the true 
solution. It has been shown that multiple points in unison may also have a stronger effect 
than an individual data point [18]. Given this, exploring the use of outlier subset 


rejection in the context of hyperbolic localization may lead to more beneficial results. 


Both source localization and collection operations presented here assume a line- 
of-sight signal path between each sensor and the transmitter. Further research is needed 
to adapt the scheme for use in a non-line-of-sight environment where there is also a 


multipath component to each signal path. 


With all wireless signal processing operations, there must be a consideration for 
the wireless communication scheme as well. Further research is required to evaluate 
signal collection from an elevated, mobile network in the context of standard digital 
wireless modulation schemes, such as binary phase-shift keying, quadrature phase-shift 
keying, and quadrature amplitude modulation. Additional effort should be devoted to the 


orthogonal frequency-division multiple access signals [85]. 
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APPENDIX. SELECTED SIMULATION CODE 


Script Name: Two-Stage Localization Versus RDOA Noise 
clc; close all; clear allj;warning off; % Program 
initialization. 


L = 10000; & Number of 
ensemble runs. 

c = 3e8; Ssignal propagation speed; 

LL = 400; Sdeployment diameter 

r = 2500; Srange from origin to source 

theta_tx = degtorad(60); %SThis is from the x-axis 

uo = [r*cos(theta_tx) r*sin(theta_tx)]'; 


% True source position. 


NN=5; %SNumber of nodes per stacks, total sensor = 3*NN 
x=[ O*ones(1,NN) -(LL/2)*ones(1,NN) (LL/2) *ones(1,NN) ]; 


je} 


%& True sensor position matrix. 
y=zeros (1,length (x) ); 
S=[x; yl; 


M = size(S,2); % Number of sensors. 
LM=length (x) ; 
N = size(S,1); 


ae 


Dimension of 


localization. 
MLE runs = 1; 
ro = sqrt (sum((uo*ones (1,M)—-S).%2))'; % True source- 
sensor ranges 
rdo = ro(2:end)-ro(l1); 
R = (eye (M-1)+ones (M-1))/2; % covariance 
structure of TDOA 
R2 = (eye (LM-1)+ones (LM-1))/2; 
NsePwrVecdB=-10:1:5; 
TSweighted least-squares= zeros(2,1); 
Linear_MLE = zeros(2,1); 
fprintf('Simulation in progress'); 
for nseIdx=1:length (NsePwrVecdB) , S vary measurement 
noise level 
Eprinte (".."); 


nsePwr = 10% (NsePwrVecdB (nseIdx) /10) 
Q = nsePwr * R; 
of TDOA (range difference) noise 
Q2 = nsePwr * R2; 
crlb (nseIdx) =sqrt (trace (TDOALOcCRLB_lin(S,uo,Q))); 
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oo Ne 


Covariance matrix 


Linear_range = 
(((1/nsePwr) *e*2*(1/ (c* (x) *%2)) %2* (LL/2)*4 -* (LM/8)) *(-1)) 57 

crlb_LM (nseIdx) =sqrt (Linear_range) ; 

Linear_bearing = (((2/nsePwr) *c*2* (1/ (c%*2) ) * (LL/2) *2 
*(LM))*{-1)); 

crlb_LM _ bearing (nseIdx)=sqrt (Linear_bearing) ; 


SimulationMSE1 = 0; 
SimulationMSE2 Or 
SimulationMSE3 = 0; 


| 
fo) 
se 


SimulationMSE1 bearing = 
SimulationMSE2_ bearing 
SimulationMSE3_bearing = 0; 


ll 
io 
se 





for k= 1 : L, % Monte Carlo 

Simulation. 

rdNse = sqrt(nsePwr/2) * randn(M,1); 

rdNsel = sqrt (nsePwr/2); 

rd = rdo + rdNse(2:end)-rdNse(1); % Noisy source 
TDOAS (range differences). 

ul = TDOALOc_LINEAR(S,rd,Q); 
[u2, senor_pos] = 
Linear_second_stage_polar_RDOA (atan2 (ul (2),ul1(1)),LM,LL/2,M 





LE_runs, [norm(uo,2),atan2 (uo(2),uo0(1))]',rdNsel, [norm(ul, 2) 
y,atan2 (ul (2),ul(1))]"); 
u2 = [u2(1)*cos(u2(2));u2(1)*sin(u2(2))]; 





[u3, senor_posl] = 
Linear_second_stage_polar_RDOA (atan2 (ul (2),ul1(1)),LM,LL/2,M 


LE_runs+3, [norm(uo,2),atan2 (uo(2),u0(1))]',rdNsel, [norm(ul, 
2), etan2 (ul ie), 1): 1?) 
u3 = [u3(1)*cos(u3(2));u3(1)*sin(u3(2))]; 
SimulationMSE1 = SimulationMSE1 + norm(ul-uo,2)%2; 
SimulationMSE2 = SimulationMSE2 + norm(u2-uo,2)%*2; 
SimulationMSE3 = SimulationMSE3 + norm(u3-uo,2)%*2; 


SimulationMSE1 bearing = SimulationMSE1_ bearing + 
(atan2 (ul (2) ,ul1(1))-theta_tx) *2; 

SimulationMSE2 bearing = SimulationMSE2_ bearing + 
(atan2 (u2(2),u2(1))-theta_tx) *2; 
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( 


S 


SimulationMSE3_bearing = SimulationMSE3_bearing + 
atan2 (u3(2),u3(1))-theta_tx) %*2; 


end; 

msel (nseIdx) = sqrt (SimulationMSE1/L) ; 

mse2 (nseIdx) = sqrt (SimulationMSE2/L) ; 

mse3 (nseIdx) = sqrt (SimulationMSE3/L) ; 
msel_bearing(nseIdx) = sqrt (SimulationMSE1_bearing/L); 


sqrt (SimulationMSE2_bearing/L); 
sqrt (SimulationMSE3_bearing/L) ; 


mse2_bearing (nselIdx) 
mse3_bearing (nselIdx) 
nd: 


fprintf('\n'); 


i 
h 


y 
1 


I 
P 


£ 


x QS O'S O'S 


yi 
1 


P 


igure(1); plot (NSePwrVecdB, (msel), 'xk', 'MarkerSize',8); 
old on; 

lot (NSsePwrVecdB, (mse2),'*r', 'MarkerSize',8); hold on; 
lot (NSePwrVecdB, (crlb),'k'); 

old on; 

lot (NsePwrVecdB, (crlb_LM),'--g'); 

rid on; hold off; 

label ('RDOA Noise STD dB(c\sigma)'); 

label ('RMSE Position (m)'); 

egend('lst Phase weighted least-squares Estimate', '2nd 
hase ML Estimate','1lst Phase CRLB RMSE Position', '2nd 








hase CRLB RMSE Position"); 

igure (3); 

lot (NSePwrVecdB, (msel_bearing),'xk', 'MarkerSize',8); hold 
nN; 

lot (NsePwrVecdB, (mse2_bearing),'*r', 'MarkerSize',8); hold 
nN; 

lot (NSePwrVecdB, (crlb_LM_ bearing), '-g'); 

rid on; hold off; 

label ('RDOA Noise STD dB(c\sigma)'); 

label('RMSE Bearing (rad)'); 








egend('lst Phase weighted least-squares Estimate', '2nd 
hase ML Estimate','2nd Phase CRLB RMSE Bearing'); 
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Script Name: Signal estimator versus wavelength for uniform position errors 


cle 

clear all 

close all 
num_run = 1000; 
rand('seed', 0) 
randn('seed', 0) 
c = 3e8; 

£ = 50:100:2600; 
f =f * le6; 
lambda = c./(f); 
k = 2*pi./lambda; 


Number_of_real_nodes_Center = 10; 
Number_of_real_nodes_Outer = Number_of_real_nodes_Center; 
ML = Number_of_real_nodes_Center*3; 

SSignal 1 
Vs_l_real = 1; 
Phase_1l_real = degtorad(0); 
DOA_1 real = 0; 


deploy_radius = 200; 





X_center = 0; 

Y_ center = 0; 

X_outer = deploy_radius; 

Y_outer = 0; 

X = 0; 

Y = 0; 

SNR = 5; 

Sscreating random receive matrix 

a = -1/10; 

b = -a; 

aa = -.9e-3; 

bb = -aa; 

Signal_estimates = zeros(6,1); 
Signal_estimates_from_mean = zeros(2,1); 
deltaX = b; 

deltaY = b; 

DOAn = bb; 

Tans_weighted least-squaresE = zeros(2,1); 


I_Q data = zeros(3,1); 
DOA_ESTIMATES = zeros(num_run,1); 
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RMSE = zeros (4, length (Number_of_real_nodes_Center) ); 
num_pts = Number_of_real_nodes_Center; 
for MCR = 1: length (k) 





MCR 

%& this is the -15 dB run 

Tans_weighted least-—squaresE = 0; 

for xx = 1l:num_run 
XnC = (a + (b-a).*rand(num_pts,1)); 
YnC = (a + (b-a).*rand(num_pts,1)); 
XnOL = (a + (b-a).*rand(num_pts,1)); 
YnOL = (a + (b-a).*rand(num_pts,1)); 
XnOR = (a + (b-a).*rand(num_pts,1)); 
YnOR = (a + (b-a).*rand(num_pts,1)); 


DOA_1l_est = Locator_Function ([2000; degtorad(70)],- 
10) - degtorad(70); 


Rx_Sig_center = Vs_1l_real.*exp((1i).* (Phase_1l_real 
+ k(MCR) .*((X_center + XnC).*sin(DOA_1l_est) + (Y_center + 
YnC) .*cos(DOA_1_est)))); 

Rx_Sig_center_w_AWGN = awgn(Rx_sig_center, SNR); 

SSummend Center Signal 

Center_IQ Est = sum(Rx_sig_center) ; 


Rx_Sig_outer_left = 
Vs_l_real.*exp((1i).* (Phase_1l_real + k(MCR) .*((- 
deploy_radius+ XnOL).*sin(DOA_1l_est) + (Y_outer + 
YnOL) .*cos(DOA_1_est)))); 

Rx_Sig_outer_left_w_AWGN = awgn(Rx_sig_outer_left, 
SNR) ; 

SEstimate Mag and Phase for Left Cluster 

Left_IQ Est = sum(Rx_sSig_outer_left); 


Rx_Sig_outer_right = 
Vs_l_real.*exp((1li).* (Phase_l_real + 
k (MCR) .* ((deploy_radius+ XnOR).*sin(DOA_l_est) + (Y_outer + 
YnOR) .*cos (DOA_l_est)))); 
Rx_Sig_outer_right_w_AWGN = 
awgn(Rx_Sig_outer_right, SNR); 
SEstimate Mag and Phase for Right Cluster 
Right_IQ_ Est = sum(Rx_sig_outer_right) ; 
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Sig_IQ_est = 
abs (sum([Center_IQ_ Est; Left_IQ_Est;Right_IQ_Est])) /ML; 








Tans_weighted least-squares = [Tans_weighted least 
squaresE Sig_IQ_ est]; 














end 
Tans_weighted least-squares = Tans_weighted least 
squaresE (2:end); 
RMSE (1,MCR) = mean(Tans_weighted least-—squaresE) ; 
Tans_weighted least-squares = 0; 
for xx = 1l:num_run 
XnC = (a + (b-a).*rand(num_pts,1)); 
YnC = (a + (b-a).*rand(num_pts,1)); 
XnOL = (a + (b-a).*rand(num_pts,1)); 
YnOL = (a + (b-a).*rand(num_pts,1)); 
XnOR = (a + (b-a).*rand(num_pts,1)); 
YnOR = (a + (b-a).*rand(num_pts,1)); 


DOA_1l_est = Locator_Function([2000; degtorad(70)],- 
5) -— degtorad(70); 


Rx_Sig_center = Vs_1l_real.*exp((1i).* (Phase_1l_real 
+ k(MCR) .*((X_center + XnC).*sin(DOA_1l_est) + (Y_center + 
YnC) .*cos(DOA_1_est)))); 

Rx_Sig_center_w_AWGN = awgn(Rx_sig_center, SNR); 

SSummend Center Signal 

Center_IQ Est = sum(Rx_sig_center) ; 


Rx_Sig_outer_left = 
Vs_l_real.*exp((1i).* (Phase_1l_real + k(MCR) .*((- 
deploy_radius+ XnOL).*sin(DOA_1l_est) + (Y_outer + 
YnOL) .*cos(DOA_1_est)))); 

Rx_Sig_outer_left_w_AWGN = awgn(Rx_sig_outer_left, 
SNR) ; 

SEstimate Mag and Phase for Left Cluster 

Left_IQ Est = sum(Rx_sig_outer_left); 


Rx_Sig_outer_right = 
Vs_l_real.*exp((1li).*(Phase_l_real + 
k (MCR) .* ((deploy_radius+ XnOR).*sin(DOA_1l_est) + (Y_outer + 
YnOR) .*cos(DOA_1_est)))); 

Rx_Sig_outer_right_w_AWGN = 
awgn (Rx_Sig_outer_right, SNR); 
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SEstimate Mag and Phase for Right Cluster 
Right_IQ Est = sum(Rx_sig_outer_right) ; 


Sig_IQ_est = 


abs (sum([Center_IQ_ Est; Left_IQ_Est;Right_IQ_Est]))/ML; 





Tans_weighted least-squaresE = [Tans_weighted 


least-squaresE Sig_IQ_ est]; 











end 

Tans_weighted least-squaresE = Tans_weighted least 
squaresE (2:end); 

RMSE (2,MCR) = mean(Tans_weighted least-—squaresE) ; 
end 

figure (1) 


plot (lambda, RMSE (1,:),lambda,RMSE(2,:)) 

xlabel ('Lambda(\lambda)'); 

ylabel('Average Normalized Gain'); 

legend('RDOA noise STD (-10 dB)','"RDOA noise STD 
dB)'); 

figure (2) 

plot (f£./le6,RMSE(1,:),£./1e6,RMSE(2,:) ) 

xlabel ('Lambda(\lambda)'); 

ylabel('Average Normalized Gain'); 

legend('RDOA noise STD (-10 dB)','RDOA noise STD 
dB)'); 
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